Package nMOLDYN :: Package Core :: Module IOFiles |
[hide private]
[frames] | no frames]

# Source Code for Module nMOLDYN.Core.IOFiles

       1  """This module implements IO-related classes, functions and procedures.
       2
       3  Classes:
       4      * TemporaryFile            : creates a temporary file stroring the evolution of an analysis.
       5      * EndOfFile                : an empty dummy class used by |DCDReader|.
       6      * FortranBinaryFile        : sets up a binary file reader.
       7      * DCDFile                  : sets up a DCD file reader.
       8      * AmberNetCDFConverter     : converts a trajectory from Amber > 9 to a MMTK NetCDF trajectory.
       9      * CHARMMConverter          : converts a trajectory from CHARMM to a MMTK NetCDF trajectory.
      10      * DL_POLYConverter         : converts a trajectory from DL_POLY > 9 to a MMTK NetCDF trajectory.
      11      * MaterialsStudioConverter : converts a trajectory from MaterialsStudio > 9 to a MMTK NetCDF trajectory.
      12      * NAMDConverter            : converts a trajectory from NAMD to a MMTK NetCDF trajectory.
      13      * VASPConverter            : converts a trajectory from VASP > 9 to a MMTK NetCDF trajectory.
      14
      15  Procedures:
      16      * convertNetCDFToASCII: converts a NetCDF file into an ASCII file.
      17      * convertASCIIToNetCDF: converts an ASCII file into a NetCDF file.
      18  """
      19
      20  # The python distribution modules
      21  from copy import copy
      22  import os
      23  import re
      24  import struct
      25  from struct import calcsize
      26  import subprocess
      27  import sys
      28  from tempfile import mktemp
      29  from time import asctime, ctime, localtime, strftime, time
      30  import xml.dom.minidom
      31
      32  # The ScientificPython modules
      33  from Scientific import N
      34  from Scientific.IO.FortranFormat import FortranFormat, FortranLine
      35  from Scientific.Geometry import Vector
      36  from Scientific.IO.NetCDF import _NetCDFFile, NetCDFFile
      37  from Scientific.IO.TextFile import TextFile
      38
      39  # The MMTK distribution modules
      40  from MMTK import Atom, AtomCluster
      41  from MMTK import Units
      42  from MMTK.ParticleProperties import Configuration, ParticleVector
      43  from MMTK.PDB import PDBConfiguration
      44  from MMTK.Trajectory import Trajectory, SnapshotGenerator, TrajectoryOutput
      45  from MMTK.Universe import InfiniteUniverse, OrthorhombicPeriodicUniverse, ParallelepipedicPeriodicUniverse
      46
      47  # The nMOLDYN modules
      48  from nMOLDYN.Core.Error import Error
      49  from nMOLDYN.Core.Mathematics import basisVectors
      50  from nMOLDYN.Core.Logger import LogMessage
      51  from nMOLDYN.Core.Preferences import PREFERENCES
      52
      53 -class TemporaryFile:
      54      """Creates a temporary file used to monitor (progress, start, end  ...) an analysis.
      55      """
      56
      57 -    def __init__(self, module = '', statusBar = None):
      58          """The constructor.
      59
      60          @param module: the name of the analysis the temporary file will be attached to.
      61          @type module: string.
      62
      63          @param statusBar: if not None, an instance of nMOLDYN.GUI.Widgets.StatusBar. The status bar, attached to the
      64              analysis, to update.
      65          @type statusBar: instance of nMOLDYN.GUI.Widgets.StatusBar
      66          """
      67
      68          self.module = module
      69          self.statusBar = statusBar
      70
      71          # retrieves the owner of the job
      72          try:
      73              import win32api
      74              owner = win32api.GetUserName()
      75
      76          except ImportError:
      77              from pwd import getpwuid
      78              from os import getuid
      79              owner = getpwuid(getuid())[0]
      80
      81          # creationTime a string containing the current time
      82          creationTime = ctime(time())
      83
      84          # pid = the pid number for the job
      85          pid = os.getpid()
      86
      87          # The suffix for the temporary file
      88          suffix = '.'.join(['',owner, str(pid), self.module, 'moldyn'])
      89
      90          # The module tempfile is used to create and handle temporary files
      91          try:
      92              from tempfile import mkstemp
      93              # This creates a unique temporary file of name 'filename' and fid 'handle'
      94              handle, self.filename = mkstemp(suffix, text=True)
      95              self.file = os.fdopen(handle, 'w')
      96
      97          except ImportError:
      98              # This just creates a unique temporary file name 'filename'
      99              self.filename = mktemp(suffix)
     100              # This open a TextFile with name 'filename'
     101              self.file     = TextFile(self.filename, 'w')
     102
     103          self.file.write(str(pid) + '\n')
     104
     105          self.file.write(creationTime + '\n')
     106
     107          self.file.write(str(owner) + '\n')
     108
     109          self.file.write(str(self.module) + '\n')
     110
     111          # Write the current time in the logfile
     112          self.file.write('0\n')
     113
     114          # Flush the internal buffer into file
     115          self.file.flush()
     116
     117          self.counter = 0
     118
     119          try:
     120              self.displayProgressRate = int(PREFERENCES.progress_rate)
     121          except ValueError:
     122              self.displayProgressRate = 10
     123
     124          if (self.displayProgressRate <= 1) or (self.displayProgressRate >= 100):
     125              self.displayProgressRate = 10
     126
     127 -    def update(self, norm):
     128          """Updates the temporary file by writing the percentage of the job that has been done.
     129
     130          @param norm: the number of steps of the outer loop of the analysis to monitor.
     131          @type norm: integer.
     132          """
     133
     134          # computes the old percentage
     135          t = int(100*self.counter/norm)
     136          oldProgress = (t/self.displayProgressRate)*self.displayProgressRate
     137          # A new step has been done
     138          self.counter += 1
     139          # computes the new percentage
     140          t = int(100*self.counter/norm)
     141          newProgress = (t/self.displayProgressRate)*self.displayProgressRate
     142
     143          if newProgress != oldProgress:
     144              LogMessage('info','Progress rate = %3d' % newProgress,['file','console'])
     145
     146          if self.statusBar is not None:
     147              self.statusBar.setValue(t)
     148
     149          # writes that percentage in the logfile. This will be useful when using the
     150          # getProcess() command.
     151          self.file.write(str(t)+"\n")
     152          # Flushes the internal buffer into file
     153          self.file.flush()
     154
     155 -    def close(self):
     156          """Closes and removes the temporary file.
     157          """
     158
     159          self.file.close()
     160          os.unlink(self.filename)
     161
     162 -class EndOfFile(Exception):
     163      pass
     164
     165 -class FortranBinaryFile(object):
     166      """Sets up a Fortran binary file reader.
     167
     168      Comments:
     169          -written by Konrad Hinsen in the scope of a DCD file reader.
     170      """
     171 -    def __init__(self, filename, byte_order = '='):
     172          self.file = file(filename, 'rb')
     173          self.byte_order = byte_order
     174
     175 -    def __iter__(self):
     176          return self
     177
     178 -    def next(self):
     179          data = self.file.read(4)
     180          if not data:
     181              raise StopIteration
     182          reclen = struct.unpack(self.byte_order + 'i', data)[0]
     183          data = self.file.read(reclen)
     184          reclen2 = struct.unpack(self.byte_order + 'i', self.file.read(4))[0]
     185          assert reclen==reclen2
     186          return data
     187
     188 -    def skipRecord(self):
     189          data = self.file.read(4)
     190          reclen = struct.unpack(self.byte_order + 'i', data)[0]
     191          self.file.seek(reclen, 1)
     192          reclen2 = struct.unpack(self.byte_order + 'i', self.file.read(4))[0]
     193          assert reclen==reclen2
     194
     195 -    def getRecord(self, format, repeat=False):
     196          try:
     197              data = self.next()
     198          except StopIteration:
     199              raise EndOfFile()
     200          if repeat:
     201              unit = struct.calcsize(self.byte_order + format)
     202              assert len(data) % unit == 0
     203              format = (len(data)/unit) * format
     204          try:
     205              return struct.unpack(self.byte_order + format, data)
     206          except:
     207              raise
     208
     209 -class DCDFile(object):
     210      """Sets up a DCD file reader.
     211      """
     212
     213 -    def __init__(self, dcd_filename):
     214          """The constructor.
     215
     216          @param dcd_filename: the name of the DCD file to read.
     217          @type dcd_filename: string.
     218          """
     219          # Identity the byte order of the file by trial-and-error
     220          self.byte_order = None
     221          data = file(dcd_filename, 'rb').read(4)
     222          for byte_order in ['<', '>']:
     223              reclen = struct.unpack(byte_order + 'i', data)[0]
     224              if reclen == 84:
     225                  self.byte_order = byte_order
     226                  break
     227          if self.byte_order is None:
     228              raise IOError("%s is not a DCD file" % dcd_filename)
     229          # Open the file
     230          self.binary = FortranBinaryFile(dcd_filename, self.byte_order)
     231          # Read the header information
     232          header_data = self.binary.next()
     233          if header_data[:4] != 'CORD':
     234              raise IOError("%s is not a DCD file" % dcd_filename)
     235          self.header = struct.unpack(self.byte_order + '9id9i', header_data[4:])
     236          self.nset = self.header[0]
     237          self.istart = self.header[1]
     238          self.nsavc = self.header[2]
     239          self.namnf = self.header[8]
     240          self.charmm_version = self.header[-1]
     241          self.has_pbc_data = False
     242          self.has_4d = False
     243          if self.charmm_version != 0:
     244              self.header = struct.unpack(self.byte_order + '9if10i',
     245                                          header_data[4:])
     246              if self.header[10] != 0:
     247                  self.has_pbc_data = True
     248              if self.header[11] != 0:
     249                  self.has_4d = True
     250          self.delta = self.header[9]*Units.akma_time
     251          # Read the title
     252          title_data = self.binary.next()
     253          nlines = struct.unpack(self.byte_order + 'i', title_data[:4])[0]
     254          assert len(title_data) == 80*nlines+4
     255          title_data = title_data[4:]
     256          title = []
     257          for i in range(nlines):
     258              title.append(title_data[:80].rstrip())
     259              title_data = title_data[80:]
     260          self.title = '\n'.join(title)
     261          # Read the number of atoms.
     262          self.natoms = self.binary.getRecord('i')[0]
     263          # Stop if there are fixed atoms.
     264          if self.namnf > 0:
     265              raise Error('NAMD converter can not handle fixed atoms yet')
     266
     267 -    def readStep(self):
     268          """Reads a frame of the DCD file.
     269          """
     270          if self.has_pbc_data:
     271              unit_cell = N.array(self.binary.getRecord('6d'), N.Float)
     272              a, gamma, b, beta, alpha, c = unit_cell
     273              if -1. < alpha < 1. and -1. < beta < 1. and -1. < gamma < 1.:
     274                  # assume the angles are stored as cosines
     275                  # (CHARMM, NAMD > 2.5)
     276                  alpha = 0.5*N.pi-N.arcsin(alpha)
     277                  beta = 0.5*N.pi-N.arcsin(beta)
     278                  gamma = 0.5*N.pi-N.arcsin(gamma)
     279              else:
     280                  # assume the angles are stored in degrees (NAMD <= 2.5)
     281                  alpha *= Units.deg
     282                  beta *= Units.deg
     283                  gamma *= Units.deg
     284              unit_cell = (a*Units.Ang, b*Units.Ang, c*Units.Ang,
     285                           alpha, beta, gamma)
     286          else:
     287              unit_cell = None
     288          format = '%df' % self.natoms
     289          x = N.array(self.binary.getRecord(format), N.Float32)*Units.Ang
     290          y = N.array(self.binary.getRecord(format), N.Float32)*Units.Ang
     291          z = N.array(self.binary.getRecord(format), N.Float32)*Units.Ang
     292          if self.has_4d:
     293              self.binary.skipRecord()
     294          return unit_cell, x, y, z
     295
     296 -    def skipStep(self):
     297          """Skips a frame of the DCD file.
     298          """
     299          nrecords = 3
     300          if self.has_pbc_data:
     301              nrecords += 1
     302          if self.has_4d:
     303              nrecords += 1
     304          for i in range(nrecords):
     305              self.binary.skipRecord()
     306
     307 -    def __iter__(self):
     308          return self
     309
     310 -    def next(self):
     311          try:
     312              return self.readStep()
     313          except EndOfFile:
     314              raise StopIteration
     315
     316 -def convertNetCDFToASCII(inputFile, outputFile, variables, floatPrecision = 9, doublePrecision = 17):
     317      """Converts a file in NetCDF format to a file in ASCII/CDL format using the ncdump program provided with
     318      the netcdf library.
     319
     320      @param inputFile: the name of the NetCDF input file.
     321      @type inputFile: string
     322
     323      @param outputFile: the name of the CDL output file.
     324      @type outputFile: string
     325
     326      @param variables: list of the NetCDF variables names (string) to extract from the NetCDF file for conversion.
     327      @type variables: list
     328
     329      @param floatPrecision: the precision on the float numbers.
     330      @type floatPrecision: integer
     331
     332      @param doublePrecision: the precision on the double numbers.
     333      @type doublePrecision: integer
     334      """
     335
     336      if not os.path.exists(inputFile):
     337          raise Error('The NetCDF input file name does not exist.')
     338
     339      if not os.path.exists(PREFERENCES.ncdump_path):
     340          raise Error('The path for ncdump program is not correctly set.')
     341
     342      if not isinstance(variables,(list,tuple)):
     343          raise Error('No NetCDF variables list given for conversion.')
     344
     345      if not variables:
     346          raise Error('Empty NetCDF variables list given for conversion.')
     347
     348      variablesList = ','.join(variables)
     349
     350      if not outputFile:
     351          raise Error('No ASCII output file name.')
     352
     353      try:
     354          output = open(outputFile, 'w')
     355          s = subprocess.call([PREFERENCES.ncdump_path, '-b', 'c', '-v',\
     356                               variablesList, '-p', '%d,%d' % (floatPrecision,doublePrecision),\
     357                               inputFile], stdout = output)
     358
     359          output.close()
     360
     361          if s:
     362              raise
     363
     364      except:
     365          raise Error('Conversion failed.')
     366
     367 -def convertASCIIToNetCDF(inputFile, outputFile):
     368      """Converts a file in ASCII format to a file in NetCDF format using the ncgen program provided with
     369      the netcdf library.
     370
     371      @param inputFile: the name of the NetCDF input file.
     372      @type inputFile: string
     373
     374      @param outputFile: the name of the CDL output file.
     375      @type outputFile: string
     376      """
     377
     378      if not os.path.exists(PREFERENCES.ncgen_path):
     379          return
     380
     381      if inputFile[-4:] == '.cdl':
     382          try:
     383              subprocess.call([PREFERENCES.ncgen_path, '-x', '-o', outputFile, inputFile])
     384          except:
     385              return
     386
     387      else:
     388          baseName = os.path.splitext(os.path.basename(inputFile))[0]
     389
     390          try:
     391              ascii = open(inputFile, 'r')
     392
     393          except IOError:
     394              raise Error('Can not open for reading %s ASCII file.' % inputFile)
     395
     396          else:
     397              asciiContents = ascii.readlines()
     398              ascii.close()
     399
     400              # Remove the empty lines of the ASCII file.
     401              asciiContents = [d.strip() for d in asciiContents if d.strip() != '']
     402
     403              globalAttributes = {}
     404              variables = []
     405
     406              globalAttributes['comment'] = ''
     407
     408              numLines = []
     409              # Loop over the non-empty lines of the ascii file.
     410              for line in asciiContents:
     411                  # If the line starts with a #, this is a comment append it to the self.comment attributes.
     412                  if line[0] == '#':
     413                      # This is the format for the declaration of a NetCDF global variable.
     414                      g = re.findall('#\s*global(.*=.*)',line.lower())
     415                      if g:
     416                          try:
     417                              gName, gVar = g[0].split('=')
     418                              globalAttributes[gName.strip()] = gVar.strip()
     419                              continue
     420
     421                          except:
     422                              raise Error('The line %s can not be parsed.' % line)
     423
     424                      v = re.findall('#\s*variable(.*)',line.lower())
     425                      if v:
     426                          try:
     427                              temp = [v.split('=') for v in v[0].strip().split(';') if v]
     428                              temp = [[r[col].strip() for r in temp] for col in range(len(temp[0]))]
     429                              d = dict(zip(temp[0],temp[1]))
     430
     431                              if not d.has_key('name'):
     432                                  raise
     433
     434                              variables.append(d)
     435                              continue
     436
     437                          except:
     438                              raise Error('The line %s can not be parsed.' % line)
     439
     440                      globalAttributes['comment'] += line[1:] + '\n'
     441
     442                  # Otherwise stores the line in its numerical form in numLines list.
     443                  else:
     444                      numLines.append([v for v in line.split()])
     445
     446              data = [[r[col].strip() for r in numLines] for col in range(len(numLines[0]))]
     447
     448              if len(variables) != len(data):
     449                  variables = [{'name' : 'column'+str(i)} for i in range(1, len(data)+1)]
     450
     451              # retrieves the owner of the job
     452              try:
     453                  import win32api
     454                  owner = win32api.GetUserName()
     455
     456              except ImportError:
     457                  from pwd import getpwuid
     458                  from os import getuid
     459                  owner = getpwuid(getuid())[0]
     460
     461              # pid = the pid number for the job
     462              pid = os.getpid()
     463
     464              suffix = '.'.join(['',owner, str(pid), 'moldyn'])
     465              # This just creates a unique temporary file name 'filename'
     466              cdlFile = mktemp(suffix)
     467
     468              file  = open(cdlFile, 'w')
     469
     470              file.write('netcdf %s {\n' % baseName)
     471
     472              file.write('dimensions:\n')
     473              file.write('\tnvalues = %d ;\n' % len(data[0]))
     474
     475              file.write('\nvariables:\n')
     476              for var in variables:
     477                  file.write('\tdouble %s(nvalues) ;\n' % var['name'])
     478                  for k, v in var.items():
     479                      if k == 'name':
     480                          continue
     481                      file.write('\t%s:%s = "%s";\n' % (var['name'],k,v))
     482
     483              file.write('\n// global attributes:\n')
     484              for k, v in globalAttributes.items():
     485                  file.write('\t:%s = "%s";\n' % (k,v))
     486
     487              file.write('\ndata:\n\n')
     488
     489              for i in range(len(variables)):
     490                  file.write(' %s =\n    ' % variables[i]['name'])
     491                  toFlush = '    '
     492                  for val in data[i]:
     493                      if len(toFlush + val + ', ') > 80:
     494                          file.write(toFlush.strip() + '\n')
     495                          toFlush = '    ' + val + ', '
     496                      else:
     497                          toFlush += val + ', '
     498
     499                  toFlush = toFlush.strip()
     500                  if toFlush[-1] == ',':
     501                      toFlush = toFlush[:-1]
     502
     503                  file.write(toFlush + ' ;\n\n')
     504
     505              file.write('}')
     506
     507              file.close()
     508
     509              try:
     510                  subprocess.call([PREFERENCES.ncgen_path, '-x', '-o', outputFile, cdlFile])
     511
     512              except:
     513                  return
     514
     515              finally:
     516                  os.unlink(cdlFile)
     517
     518 -class AmberNetCDFConverter:
     519      """Converts an Amber NetCDF Trajectory into a MMTK NetCDFFile.
     520
     521      Comments:
     522          - this code is an improved version of the original converter written by Paolo Calligari.
     523      """
     524
     525 -    def __init__(self, pdbFile, amberNetCDFFile, outputFile, timeStep = 1.0):
     526          """
     527          The constructor. Will do the conversion.
     528
     529          @param pdbFile: the Amber PDB file name of a frame of the trajectory to convert.
     530          @type pdbFile: string
     531
     532          @param amberNetCDFFile: the Amber NetCDF file name of the trajectory to convert.
     533          @type amberNetCDFFile: string
     534
     535          @param outputFile: the name of MMTK NetCDF trajectory output file.
     536          @type outputFile: string
     537
     538          @param timeStep: the timestep that will used when building the MMTK trajectory. Default to 1 ps.
     539          @type timeStep: float.
     540          """
     541
     542          self.pdbFile = pdbFile
     543          self.amberNetCDFFile = amberNetCDFFile
     544          self.outputFile = outputFile
     545          self.timeStep = timeStep
     546
     547          # Do the conversion.
     548          self.__convert()
     549
     550 -    def __convert(self):
     551          """Performs the actual conversion.
     552          """
     553
     554          # The Amber NetCDF input file is opened for reading.
     555          try:
     556              netcdfInputFile = NetCDFFile(self.amberNetCDFFile, 'r')
     557
     558          except:
     559              raise Error('The Amber NetCDF file %s could not be opened for reading.' % self.amberNetCDFFile)
     560
     561          else:
     562
     563              nSteps = len(netcdfInputFile.variables['time'])
     564
     565              # Case of a periodic universe.
     566              if netcdfInputFile.variables.has_key('cell_lengths'):
     567                  # The box dimensions are stored in |box| and converted into nanometers.
     568                  box = netcdfInputFile.variables['cell_lengths']
     569
     570                  # The box angles are stored in |angles| and converted into radians.
     571                  angles = netcdfInputFile.variables['cell_angles']
     572
     573                  cellParams = list(box[0]*Units.Ang) + list(angles[0]*Units.deg)
     574
     575                  universe = ParallelepipedicPeriodicUniverse(basisVectors(cellParams))
     576
     577              else:
     578                  universe = InfiniteUniverse()
     579
     580          # The PDB input file is opened for reading.
     581          try:
     582              # Get data about the universe.
     583              conf = PDBConfiguration(self.pdbFile)
     584
     585          except:
     586              raise Error('The PDB file %s could not be opened for reading.' % self.pdbFile)
     587
     588          else:
     589              molecules = conf.createAll()
     590              universe.addObject(molecules)
     591
     592          trajectory = Trajectory(universe, self.outputFile, 'w', 'Converted from AMBER NetCDF file.')
     593          snapshot = SnapshotGenerator(universe, actions=[TrajectoryOutput(trajectory, ['all'], 0, None, 1)])
     594
     595          num = molecules.numberOfAtoms()
     596
     597          for i in range(nSteps):
     598
     599              t = float(i+1)*self.timeStep
     600
     601              # Case of a periodic universe. The cell parameters are calculated.
     602              if universe.is_periodic:
     603                  cellParams = list(box[i]*Units.Ang) + list(angles[i]*Units.deg)
     604                  b = basisVectors(cellParams)
     605                  cell = N.zeros((9,), typecode = N.Float)
     606                  cell[0:3] = b[0]
     607                  cell[3:6] = b[1]
     608                  cell[6:9] = b[2]
     609
     610              else:
     611                  cell = None
     612
     613              # The coordinates are converted in nanometers.
     614              coordinates = netcdfInputFile.variables['coordinates'][i][:num]*Units.Ang
     615
     616              # The universe configuration is updated.
     617              universe.setConfiguration(Configuration(universe, coordinates, cell))
     618
     619              # The frame is written.
     620              snapshot(data = {'time': t})
     621
     622          trajectory.close()
     623          netcdfInputFile.close()
     624
     625 -class CHARMMConverter:
     626      """Converts a CHARMM Trajectory into a MMTK NetCDFFile.
     627
     628      Comments:
     629          - this code is based on the original converter written by Konrad Hinsen.
     630      """
     631
     632 -    def __init__(self, pdbFile, dcdFile, outputFile):
     633          """
     634          The constructor. Will do the conversion.
     635
     636          @param pdbFile: the CHARMM PDB file name of a frame of the trajectory to convert.
     637          @type pdbFile: string
     638
     639          @param dcdFile: the CHARMM DCD file name of the trajectory to convert.
     640          @type dcdFile: string
     641
     642          @param outputFile: the name of MMTK NetCDF trajectory output file.
     643          @type outputFile: string
     644          """
     645
     646          self.pdbFile = pdbFile
     647          self.dcdFile = dcdFile
     648          self.outputFile = outputFile
     649
     650          # Do the conversion.
     651          self.__convert()
     652
     653 -    def __convert(self):
     654          """Performs the actual conversion.
     655          """
     656
     657          # Open the DCD trajectory file
     658          dcd = DCDFile(self.dcdFile)
     659          step = dcd.istart
     660          step_increment = dcd.nsavc
     661          dt = dcd.delta
     662          unit_cell, x, y, z = dcd.readStep()
     663
     664          # The DCD file does not contain any periodic boundary condition data.
     665          if not dcd.has_pbc_data:
     666              universe = InfiniteUniverse()
     667          # The DCD file contains periodic boundary condition datas. Create a parallelepipedic universe out of
     668          # the first unit cell defined for the first step of the trajectory.
     669          else:
     670              universe = ParallelepipedicPeriodicUniverse(basisVectors(unit_cell))
     671
     672          # Create all objects from the PDB file. The PDB file must match the
     673          # the DCD file (same atom order).
     674          conf = PDBConfiguration(self.pdbFile)
     675          universe.addObject(conf.createAll())
     676
     677          # Create the output trajectory
     678          trajectory = Trajectory(universe, self.outputFile, mode = 'w', comment = dcd.title)
     679          snapshot = SnapshotGenerator(universe, actions=[TrajectoryOutput(trajectory, ["all"], 0, None, 1)])
     680
     681          # Iterate over the DCD file and write to the output trajectory
     682          conf_array = universe.configuration().array
     683          comp = 0
     684          while True:
     685              conf_array[:, 0] = x
     686              conf_array[:, 1] = y
     687              conf_array[:, 2] = z
     688
     689              if universe.is_periodic:
     690                  universe.setShape(basisVectors(unit_cell))
     691
     692              step_data = {'time': step*dt, 'step': step}
     693              snapshot(data = step_data)
     694              step += step_increment
     695
     696              try:
     697                  unit_cell, x, y, z = dcd.readStep()
     698                  comp += 1
     699
     700              except EndOfFile:
     701                  break
     702
     703          # Close the output trajectory
     704          trajectory.close()
     705
     706 -class DL_POLYConverter:
     707      """Converts a DL_POLY Trajectory into a MMTK NetCDFFile.
     708      """
     709
     710 -    def __init__(self, fieldFile, historyFile, outputFile, specialAtoms = {}):
     711          """
     712          The constructor. Will do the conversion.
     713
     714          @param fieldFile: the DL_POLY FIELD file name of the trajectory to convert.
     715          @type fieldFile: string
     716
     717          @param historyFile: the DL_POLY HISTORY file name of the trajectory to convert.
     718          @type historyFile: string
     719
     720          @param outputFile: the name of MMTK NetCDF trajectory output file.
     721          @type outputFile: string
     722
     723          @param specialAtoms: dictionnary of the form {s1 : e1, s2 : e2 ...} where 's1', 's2' ... and 'e1', 'e2' ...
     724              are respectively the DL_POLY name and the symbol of atoms 1, 2 ...
     725          @type specialAtoms: dict
     726          """
     727
     728          self.fieldFile = fieldFile
     729          self.historyFile = historyFile
     730          self.outputFile = outputFile
     731
     732          self.specialAtoms = copy(specialAtoms)
     733
     734          # All the items of the special atoms dictionnary are lowerized.
     735          for k,v in self.specialAtoms.items():
     736              if k.lower() != k:
     737                  self.specialAtoms[k.lower()] = v.lower()
     738                  del self.specialAtoms[k]
     739              else:
     740                  self.specialAtoms[k] = v.lower()
     741
     742          # Do the conversion.
     743          self.__convert()
     744
     745 -    def __readFIELD(self):
     746          """Reads the DL_POLY FIELD file that contains informations about the system.
     747          """
     748
     749          nspecies = 0
     750
     751          # The FIELD file is opened for reading and stored into |lines|.
     752          fieldFile = open(self.fieldFile, 'r')
     753          lines = fieldFile.readlines()
     754          fieldFile.close()
     755
     756          while 1:
     757              # Check whether the FIELD keyword 'molecules' or 'molecular type' have been matched.
     758              if re.match('MOLECULES', lines[0].upper()):
     759                  # |nspecies| gives the number of different type of molecules in the system.
     760                  nspecies = int(lines[0].split()[1])
     761                  break
     762
     763              if re.match('MOLECULAR TYPES', lines[0].upper()):
     764                  # |nspecies| gives the number of different type of molecules in the system.
     765                  nspecies = int(lines[0].split()[2])
     766                  break
     767              lines = lines[1:]
     768          lines = lines[1:]
     769
     770          if nspecies == 0:
     771              raise Error('No molecular type was found in your DL_POLY/FIELD file.')
     772
     773          # |species| is a tuplet list that gives the main informations about the molecular system.
     774          species = []
     775          for i in range(nspecies):
     776              # |name| is the name of the molecular type.
     777              name = lines[0].strip()
     778              # n is the number of molecules of |name| molecular type.
     779              n = int(lines[1].split()[1])
     780
     781              # |natoms| is the number of atoms in molecule |name|.
     782              natoms = int(lines[2].split()[1])
     783
     784              lines = lines[3:]
     785              atoms = []
     786              while natoms > 0:
     787                  # Each atom line is splitted in order to fetch the different atomic parameter
     788                  # Take care the second and third values must be real, while the fourth and
     789                  # fifth ones, although optional, must be integer.
     790                  data = lines[0].split()
     791
     792                  if len(data) < 3:
     793                      raise Error('The atomic record must contain at least three element.')
     794
     795                  lines = lines[1:]
     796                  atom_name = data[0].strip().lower()
     797
     798                  element = atom_name
     799                  # First of all checks whether the atom has been defined as a special atoms.
     800                  if element in self.specialAtoms.keys():
     801                      element = self.specialAtoms[element]
     802
     803                  else:
     804
     805                      element = atom_name
     806                      while len(element) > 0:
     807                          try:
     808                              Atom(element)
     809
     810                          except:
     811                              element = element[:-1]
     812
     813                          else:
     814                              break
     815
     816                      if not element:
     817                          raise Error('The element corresponding to atom %s could not be found in the MMTK Database.' % atom_name)
     818
     819                  try:
     820                      atomweight = float(data[1])
     821
     822                  except:
     823                      raise Error('The second value of the atomic record (atom weight) must be a number.')
     824
     825                  try:
     826                      atomchg = float(data[2])
     827
     828                  except:
     829                      raise Error('The third value of the atomic record (atom charge) must be a number.')
     830
     831                  try:
     832                      nrepeat = max(int(data[3]), 1)
     833
     834                  except:
     835                      nrepeat = 1
     836
     837                  for j in range(nrepeat):
     838                      atoms.append((element, atom_name))
     839
     840                  natoms = natoms - nrepeat
     841
     842              # Check for keyword 'FINISH' or 'CONSTRAINTS' to stop.
     843              while ((not re.match('FINISH', lines[0].upper())) and (not re.match('CONSTRAINTS', lines[0].upper()))):
     844                  lines = lines[1:]
     845
     846              # |constraints| is a tuplet list of the bond contraints found in the molecule.
     847              # Its element are of the form (atom 1 index, atom 2 index, bond length).
     848              constraints = []
     849              if re.match('CONSTRAINTS',lines[0].upper()):
     850                  # |nc| is the number of contrained bonds in the molecule
     851                  nc = int(lines[0].split()[1])
     852                  lines = lines[1:]
     853
     854                  while nc > 0:
     855                      l = lines[0].split()
     856                      i1 = int(l[0])-1
     857                      i2 = int(l[1])-1
     858                      d = float(l[2])*Units.Ang
     859                      constraints.append((i1, i2, d))
     860                      lines = lines[1:]
     861                      nc = nc - 1
     862
     863                  while not re.match('FINISH', lines[0].upper()):
     864                      lines = lines[1:]
     865
     866              lines = lines[1:]
     867              species.append([name, n, atoms, constraints])
     868
     869          return species
     870
     871 -    def __readHISTORY(self):
     872          """Reads the HISTORY file and write the output trajectory frame by frame.
     873          """
     874
     875          historyFile = open(self.historyFile, 'r')
     876
     877          # This is the first line of the HISTORY file that contains informations about the system.
     878          header = historyFile.readline().strip()
     879
     880          # The second line of the HISTORY file gives some infos about the trajectory.
     881          # There must be a maximum of 3 integers.
     882          # The first one refers to the trajectory key (O = only coordinates, 1 = only velocities, 2 = both)
     883          # The second one refers to the periodic boundariy key.
     884          keytrj, imcon, natms = FortranLine(historyFile.readline(), '3I10')
     885          if imcon >= 4:
     886              raise Error('Box shape not supported by MMTK')
     887
     888          # The periodic boundary key of 0 refers to non periodic boundary in DL_POLY
     889          if imcon == 0:
     890              universe = InfiniteUniverse()
     891
     892          else:
     893              universe = ParallelepipedicPeriodicUniverse()
     894
     895          for mol_name, mol_count, atoms, constraints in self.molecules:
     896              for i in range(mol_count):
     897                  number = 0
     898                  atom_objects = []
     899                  for element, name in atoms:
     900                      a = Atom(element, name = name+str(number))
     901                      a.number = number
     902                      number = number + 1
     903                      atom_objects.append(a)
     904
     905                  if len(atom_objects) == 1:
     906                      universe.addObject(atom_objects[0])
     907
     908                  else:
     909                      # When the molecular system is made of more than one atom. They are first gathered
     910                      # into an MMTK AtomCluster that is added to the universe.
     911                      ac = AtomCluster(atom_objects, name = mol_name)
     912                      for i1, i2, d in constraints:
     913                          ac.addDistanceConstraint(atom_objects[i1], atom_objects[i2], d)
     914                      universe.addObject(ac)
     915
     916          if keytrj > 0:
     917              universe.initializeVelocitiesToTemperature(0.)
     918
     919          # The trajectory is opened for writing.
     920          trajectory = Trajectory(universe, self.outputFile, mode = 'w', comment = header)
     921
     922          # Each time snapshot is called, the universe contents i flushed into the output file.
     923          snapshot = SnapshotGenerator(universe, actions = [TrajectoryOutput(trajectory, ["all"], 0, None, 1)])
     924
     925          conf = universe.configuration()
     926          vel = universe.velocities()
     927
     928          # |grad| is a ParticleVector MMTK object that associates a vector to each element of the universe.
     929          # That vector is the force associated to that element at a given time step.
     930          # This should not be of any use in nMoldyn but, as DLPOLY can provide it, it is done just in case ...
     931          grad = ParticleVector(universe)
     932
     933          # Some predined format.
     934          history_timestep_line = FortranFormat('A8,4I10,F12.6')
     935          history_pbc_line = FortranFormat('3G12.4')
     936
     937          while 1:
     938              # The HISTORY file was not closed so the reading continue from the third line.
     939              line = historyFile.readline()
     940              # An empty line means the EOF.
     941              if not line:
     942                  break
     943
     944              # The third line contains some time informations about the trajectory.
     945              data = FortranLine(line, history_timestep_line)
     946
     947              # The current time step.
     948              nstep = data[1]
     949
     950              # The number of atoms in the configuration.
     951              natoms = data[2]
     952
     953              keytrj = data[3]
     954
     955              imcon = data[4]
     956
     957              # The trajectory time step.
     958              tstep = data[5]
     959
     960              step_data = {'time': nstep*tstep}
     961
     962              if imcon:
     963                  # This line refers to the basis vector 'a' in cartesian coordinates.
     964                  data = FortranLine(historyFile.readline(), history_pbc_line)
     965                  dirA = Vector(data)*Units.Ang
     966
     967                  # This line refers to the basis vector 'b' in cartesian coordinates.
     968                  data = FortranLine(historyFile.readline(), history_pbc_line)
     969                  dirB = Vector(data)*Units.Ang
     970
     971                  # This line refers to the basis vector 'c' in cartesian coordinates.
     972                  data = FortranLine(historyFile.readline(), history_pbc_line)
     973                  dirC = Vector(data)*Units.Ang
     974
     975                  cell = N.zeros((9,), typecode = N.Float)
     976
     977                  cell[0:3] = dirA
     978                  cell[3:6] = dirB
     979                  cell[6:9] = dirC
     980
     981              else:
     982                  cell = None
     983
     984              for i in range(natoms):
     985                  historyFile.readline()
     986                  # The configuration of atom |i| is updated.
     987                  conf.array[i] = [float(v) for v in historyFile.readline().split()]
     988                  # Case where velocities are given in the HISTORY file.
     989                  if keytrj > 0:
     990                      vel.array[i] = [float(v) for v in historyFile.readline().split()]
     991                      # Case where forces are also given in the HISTORY file.
     992                      if keytrj > 1:
     993                          grad.array[i] = [float(v) for v in historyFile.readline().split()]
     994
     995              universe.setConfiguration(Configuration(universe, conf.array*Units.Ang, cell))
     996
     997              if keytrj > 0:
     998                  universe.setVelocities(vel*Units.Ang/Units.ps)
     999
    1000              if keytrj > 1:
    1001                  N.multiply(grad.array, -Units.amu*Units.Ang/Units.ps**2, grad.array)
    1002                  step_data['gradients'] = grad
    1003
    1004              snapshot(data = step_data)
    1005
    1006          trajectory.close()
    1007
    1008          historyFile.close()
    1009
    1010 -    def __convert(self):
    1011          """Performs the actual conversion.
    1012          """
    1013
    1014          # |self.molecules| is a list that contains the main informations about each molecular type found in the system.
    1015          # Each element of that list corresponds to a molecular type found in the system and is the tuplet that
    1016          # contains respectively the name of the molecular type, the number of molecules of that molecular type, the
    1017          # list of atoms of one molecule, and the bond of the molecule with a constraints.
    1018          self.molecules = self.__readFIELD()
    1019
    1020          self.__readHISTORY()
    1021
    1022 -class MaterialsStudioConverter:
    1023      """Converts a MaterialsStudio Discover or Forcite Trajectory into a MMTK NetCDFFile.
    1024      """
    1025
    1026      atomLineFormat = FortranFormat('A5,1X,F14.9,1X,F14.9,1X,F14.9,1X,A4,1X,A7,A7,1X,A2,1X,F6.3')
    1027
    1028 -    def __init__(self, module, xtdxsdFile, histrjFile, outputFile, subselection = None):
    1029          """
    1030          The constructor. Will do the conversion.
    1031
    1032          @param module: a string being one of 'Discover' or 'Forcite' specifying which module of MaterialsStudio the
    1033              trajectory is coming from.
    1034          @type module: string
    1035
    1036          @param xtdxsdFile: the MaterialsStudio XTD or XSD file name of the trajectory to convert.
    1037          @type xtdxsdFile: string
    1038
    1039          @param histrjFile: the MaterialsStudio HIS (Discover) or TRJ (Forcite) file name of the trajectory to convert.
    1040          @type histrjFile: string
    1041
    1042          @param outputFile: the name of MMTK NetCDF trajectory output file.
    1043          @type outputFile: string
    1044
    1045          @param subselection: if not None, list of the indexes (integer >= 1) of the atoms to select when writing out
    1046              the MMTK trajectory. The order being the one defined in the XTD/XSD file.
    1047          @type subselection: list
    1048          """
    1049
    1050          self.module = module
    1051          self.xtdxsdFile = xtdxsdFile
    1052          self.histrjFile = histrjFile
    1053          self.outputFile = outputFile
    1054          self.subselection = subselection
    1055
    1056          # Do the conversion.
    1057          self.__convert()
    1058
    1059 -    def createCluster(self, at, clust):
    1060
    1061          if at.has_key('cluster'):
    1062              return
    1063
    1064          at['cluster'] = True
    1065
    1066          clust.append(at)
    1067
    1068          for cAt in at['connections']:
    1069              self.createCluster(self.atomContents[self.atomIndexes[cAt]], clust)
    1070
    1071 -    def readXTDFile(self):
    1072          """Reads the Materials Studio XTD or  XSD file and set up the universe from which the NetCDF
    1073          MMTK trajectory will be written.
    1074
    1075          @note: the XTD and XSD file are xml file.
    1076          """
    1077
    1078          xmlFile = xml.dom.minidom.parse(self.xtdxsdFile)
    1079
    1080          identityMapping = xmlFile.getElementsByTagName('IdentityMapping')[0]
    1081
    1082          atomsTag = identityMapping.getElementsByTagName('Atom3d')
    1083
    1084          try:
    1085              spaceGroup = identityMapping.getElementsByTagName('SpaceGroup')[0]
    1086
    1087              self.normalizedCell = N.zeros((9,), typecode = N.Float)
    1088              self.normalizedCell[0:3] = Vector([float(v) for v in spaceGroup.getAttribute('AVector').split(',')]).normal()
    1089              self.normalizedCell[3:6] = Vector([float(v) for v in spaceGroup.getAttribute('BVector').split(',')]).normal()
    1090              self.normalizedCell[6:9] = Vector([float(v) for v in spaceGroup.getAttribute('CVector').split(',')]).normal()
    1091
    1092              self.universe = ParallelepipedicPeriodicUniverse()
    1093
    1094          except:
    1095              LogMessage('No SpaceGroup tag found in the xtd/xsd file. The MMTK universe created will be an infinite universe.')
    1096              self.universe = InfiniteUniverse()
    1097
    1098          self.clusters = []
    1099          if self.subselection is None:
    1100
    1101              bondsTag = xmlFile.getElementsByTagName('IdentityMapping')[0].getElementsByTagName('Bond')
    1102              bondsContents = {}
    1103              for bd in bondsTag:
    1104                  bdIndex = int(bd.getAttribute('ID'))
    1105                  bondsContents[bdIndex] = [int(v) for v in bd.getAttribute('Connects').split(',')]
    1106
    1107              self.atomIndexes = {}
    1108              self.atomContents = []
    1109              comp = 0
    1110              for at in atomsTag:
    1111                  atId = int(at.getAttribute('ID'))
    1112                  element = str(at.getAttribute('Components')).split(',')[0].strip()
    1113                  coord = N.array([float(v) for v in at.getAttribute('XYZ').split(',')])*Units.Ang
    1114
    1115                  connections = []
    1116                  if at.hasAttribute('Connections'):
    1117                      connectedBonds = [int(v) for v in at.getAttribute('Connections').split(',')]
    1118                      for b in connectedBonds:
    1119                          if bondsContents.has_key(b):
    1120                              for bb in bondsContents[b]:
    1121                                  if bb != atId:
    1122                                      connections.append(bb)
    1123
    1124                  data = {'coord' : coord, 'element' : element, 'connections' : connections, 'serial_number' : comp}
    1125                  self.atomContents.append(data)
    1126
    1127                  self.atomIndexes[atId] = comp
    1128
    1129                  comp +=1
    1130
    1131              for at in self.atomContents:
    1132
    1133                  if at.has_key('cluster'):
    1134                      continue
    1135
    1136                  else:
    1137                      clust = []
    1138                      self.createCluster(at, clust)
    1139                      self.clusters.append(clust)
    1140
    1141          else:
    1142
    1143              for cInfo in self.subselection:
    1144
    1145                  subCluster = []
    1146                  for ind in cInfo:
    1147                      atXML = atomsTag[ind-1]
    1148                      element = str(atXML.getAttribute('Components')).split(',')[0].strip()
    1149                      data = {'coord'         : N.array([float(v) for v in atXML.getAttribute('XYZ').split(',')])*Units.Ang,
    1150                              'element'       : element,
    1151                              'serial_number' : ind - 1}
    1152                      subCluster.append(data)
    1153                  self.clusters.append(subCluster)
    1154
    1155          xmlFile.unlink()
    1156
    1157          for c in range(len(self.clusters)):
    1158              clust = self.clusters[c]
    1159              atoms = []
    1160              atomNames = {}
    1161              for atom in clust:
    1162                  if atomNames.has_key(atom['element']):
    1163                      atomNames[atom['element']] += 1
    1164                  else:
    1165                      atomNames[atom['element']] = 1
    1166                  atoms.append(Atom(atom['element'], name = str(atomNames[atom['element']])))
    1167
    1168              clustName = ''
    1169              for k, v in atomNames.items():
    1170                  clustName += '%s%d' % (k,v)
    1171              clustName += '_' + str(c+1)
    1172              self.universe.addObject(AtomCluster(atoms, name = clustName))
    1173
    1174 -    def readHISFile(self):
    1175          """Reads a Materials Studio HIS file and fills up the NetCDF trajectory file.
    1176          """
    1177
    1178          selCoordinates = N.zeros((self.universe.numberOfAtoms(),3), typecode = N.Float)
    1179          selVelocities = N.zeros((self.universe.numberOfAtoms(),3), typecode = N.Float)
    1180
    1181          # The output trajectory NetCDF file is opened for writing.
    1182          trajectory = Trajectory(self.universe, self.outputFile, mode = 'w', comment = 'From Discover trajectory file')
    1183
    1184          snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(trajectory,["all"], 0, None, 1)])
    1185
    1186          # Read the his file.
    1187          hisFile = open(self.histrjFile,'rb')
    1188
    1189          # Record 1
    1190          rec = '!4xi8x'
    1191          recSize = calcsize(rec)
    1192          hisFile.read(recSize)
    1193
    1194          # Record 2
    1195          rec = '!' + '4s' * 20 + 'd8x'
    1196          recSize = calcsize(rec)
    1197          data = struct.unpack(rec, hisFile.read(recSize))
    1198          Vershn = data[20]
    1199
    1200          LogMessage('info','Vershn -- > %s' % Vershn,['file'])
    1201
    1202          # Record 3
    1203          rec = '!' + '4s' * 20 + '8x'
    1204          recSize = calcsize(rec)
    1205          hisFile.read(recSize)
    1206
    1207          # Record 4
    1208          rec = '!' + '4s' * 20 + '8x'
    1209          recSize = calcsize(rec)
    1210          hisFile.read(recSize)
    1211
    1212          # Record 5
    1213          rec = '!i'
    1214          recSize = calcsize(rec)
    1215          data = struct.unpack(rec, hisFile.read(recSize))
    1216          NAtTyp = data[0]
    1217
    1218          rec = '!' + '4s' * NAtTyp + '%sd8x' % NAtTyp
    1219          recSize = calcsize(rec)
    1220          hisFile.read(recSize)
    1221
    1222          # Record 6
    1223          rec = '!i'
    1224          recSize = calcsize(rec)
    1225          data = struct.unpack(rec, hisFile.read(recSize))
    1226          NNmRes = data[0]
    1227
    1228          rec = '!' + '4s' * NNmRes + '8x'
    1229          recSize = calcsize(rec)
    1230          hisFile.read(recSize)
    1231
    1232          # Record 7
    1233          rec = '!i'
    1234          recSize = calcsize(rec)
    1235          data = struct.unpack(rec, hisFile.read(recSize))
    1236          NAtoms = data[0]
    1237
    1238          rec = '!%si' % NAtoms
    1239          recSize = calcsize(rec)
    1240          hisFile.read(recSize)
    1241
    1242          if Vershn < 2.9:
    1243              rec = '!' + '4s' * NAtoms
    1244
    1245          else:
    1246              rec = '!' + '5s' * NAtoms
    1247
    1248          recSize = calcsize(rec)
    1249          hisFile.read(recSize)
    1250
    1251          hisFile.read(calcsize('!8x'))
    1252
    1253          # Record 8
    1254          rec = '!ii'
    1255          recSize = calcsize(rec)
    1256          data = struct.unpack(rec, hisFile.read(recSize))
    1257
    1258          NAtMov = data[1]
    1259
    1260          if Vershn >= 2.6:
    1261              rec = '!%si' % NAtMov
    1262              recSize = calcsize(rec)
    1263              data = struct.unpack(rec, hisFile.read(recSize))
    1264              movableatoms = data
    1265              LogMessage('info','movableatoms -- > %s' % str(movableatoms),['file'])
    1266
    1267          hisFile.read(calcsize('!8x'))
    1268
    1269          # Record 9
    1270          rec = '!i'
    1271          recSize = calcsize(rec)
    1272          data = struct.unpack(rec, hisFile.read(recSize))
    1273          NMol = data[0]
    1274
    1275          rec = '!%si%si' % (NMol,NMol)
    1276          recSize = calcsize(rec)
    1277          hisFile.read(recSize)
    1278
    1279          hisFile.read(calcsize('!8x'))
    1280
    1281          # Record 10
    1282          rec = '!i'
    1283          recSize = calcsize(rec)
    1284          data = struct.unpack(rec, hisFile.read(recSize))
    1285          NRes = data[0]
    1286
    1287          rec = '!%si%si' % (NRes*2,NRes)
    1288          recSize = calcsize(rec)
    1289          hisFile.read(recSize)
    1290
    1291          hisFile.read(calcsize('!8x'))
    1292
    1293          # Record 11
    1294          rec = '!i'
    1295          recSize = calcsize(rec)
    1296          data = struct.unpack(rec, hisFile.read(recSize))
    1297          NBonds = data[0]
    1298
    1299          if NBonds > 0 :
    1300              rec = '!%si' % 2*NBonds
    1301              recSize = calcsize(rec)
    1302              hisFile.read(recSize)
    1303
    1304          hisFile.read(calcsize('!8x'))
    1305
    1306          # Record 12
    1307          rec = '!4149di4di6d6i8x'
    1308          recSize = calcsize(rec)
    1309          data = struct.unpack(rec, hisFile.read(recSize))
    1310
    1311          # Record 13
    1312          rec = '!idii8x'
    1313          recSize = calcsize(rec)
    1314          data = struct.unpack(rec, hisFile.read(recSize))
    1315          NEner, timestep, frequency, startingstep = data
    1316
    1317          LogMessage('info','NEner -- > %s' % NEner,['file'])
    1318          LogMessage('info','timestep -- > %s fs' % timestep,['file'])
    1319          LogMessage('info','frequency -- > %s step(s)' % frequency,['file'])
    1320          LogMessage('info','startingstep -- > %s' % startingstep,['file'])
    1321
    1322          # Record 14
    1323          nrjBlockSize = 59+5*NMol+NEner+NMol*NEner
    1324          rec = '!%sd8x' % nrjBlockSize
    1325          recSize = calcsize(rec)
    1326          data = struct.unpack(rec, hisFile.read(recSize))
    1327          totalenergy = data[0]
    1328
    1329          LogMessage('info','totalenergy -- > %s' % totalenergy,['file'])
    1330
    1331          nAtoms3 = 3*NAtoms
    1332
    1333          # Record 15
    1334          rec = '!%sf8x' % nAtoms3
    1335          recSize = calcsize(rec)
    1336          data = struct.unpack(rec, hisFile.read(recSize))
    1337          firstatomxyz = data[0:3]
    1338
    1339          LogMessage('info','firstatomxyz -- > %s' % str(firstatomxyz),['file'])
    1340
    1341          # Record 16
    1342          rec = '!%sf8x' % nAtoms3
    1343          recSize = calcsize(rec)
    1344          data = struct.unpack(rec, hisFile.read(recSize))
    1345          firstatomvel = data[0:3]
    1346
    1347          LogMessage('info','firstatomvel -- > %s' % str(firstatomvel),['file'])
    1348
    1349          # Frame record
    1350          frame = 0
    1351          while True:
    1352
    1353              try:
    1354                  # Record N
    1355                  rec = '!i8x'
    1356                  recSize = calcsize(rec)
    1357                  hisFile.read(recSize)
    1358
    1359                  # Record N + 1
    1360                  rec = '!%sd8x' % nrjBlockSize
    1361                  recSize = calcsize(rec)
    1362                  hisFile.read(recSize)
    1363
    1364                  # Record N + 2
    1365                  rec = '!6d9d8x'
    1366                  recSize = calcsize(rec)
    1367                  data = struct.unpack(rec, hisFile.read(recSize))
    1368
    1369                  cell = N.array(data[6:])*Units.Ang
    1370                  LogMessage('info','cell -- > %s' % str(cell),['file'])
    1371
    1372                  # Record N+3
    1373                  rec = '!%sf8x' % nAtoms3
    1374                  recSize = calcsize(rec)
    1375                  data = struct.unpack(rec, hisFile.read(recSize))
    1376                  xyz = N.array(data).reshape((NAtoms,3))*Units.Ang
    1377
    1378                  # Record N+4
    1379                  rec = '!%sf' % nAtoms3
    1380                  recSize = calcsize(rec)
    1381                  data = struct.unpack(rec, hisFile.read(recSize))
    1382                  vel = N.array(data).reshape((NAtoms,3))*Units.Ang
    1383
    1384                  hisFile.read(calcsize('!8x'))
    1385
    1386              except:
    1387                  break
    1388
    1389              else:
    1390                  frame += 1
    1391
    1392                  comp = 0
    1393                  for clust in self.clusters:
    1394                      for subClust in clust:
    1395                          selCoordinates[comp,:] = xyz[subClust['serial_number'],:]
    1396                          selVelocities[comp,:] = vel[subClust['serial_number'],:]
    1397                          comp += 1
    1398
    1399                  if not self.universe.is_periodic:
    1400                      cell = None
    1401                  self.universe.setConfiguration(Configuration(self.universe, selCoordinates, cell))
    1402                  self.universe.setVelocities(ParticleVector(self.universe, selVelocities))
    1403
    1404                  snapshot(data = {'time': (startingstep + frequency*frame)*timestep*Units.fs})
    1405
    1406          # The his file is closed.
    1407          hisFile.close()
    1408
    1409          trajectory.close()
    1410
    1411 -    def readTRJFile(self):
    1412          """Reads a Materials Studio HIS file and fills up the NetCDF trajectory file.
    1413          """
    1414
    1415          selCoordinates = N.zeros((self.universe.numberOfAtoms(),3), typecode = N.Float)
    1416          selVelocities = N.zeros((self.universe.numberOfAtoms(),3), typecode = N.Float)
    1417
    1418          # The output trajectory NetCDF file is opened for writing.
    1419          trajectory = Trajectory(self.universe, self.outputFile, mode = 'w', comment = 'From Forcite trajectory file')
    1420
    1421          snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(trajectory,["all"], 0, None, 1)])
    1422
    1423          # Read the his file.
    1424          trjFile = open(self.histrjFile,'rb')
    1425
    1426          # Record 1
    1427          rec = '!4x4s20i8x'
    1428          recSize = struct.calcsize(rec)
    1429          data = struct.unpack(rec, trjFile.read(recSize))
    1430
    1431          HDR = data[0]
    1432          LogMessage('info','HDR -- > %s' % HDR,['file'])
    1433
    1434          ICNTRL = data[1:]
    1435          LogMessage('info','ICNTRL -- > %s' % str(ICNTRL),['file'])
    1436
    1437          if ICNTRL[0] == 2000:
    1438              prec = 'f'
    1439
    1440          elif ICNTRL[0] == 2010:
    1441              prec = 'd'
    1442
    1443          elif ICNTRL[0] == 3000:
    1444              prec = 'd'
    1445
    1446          else:
    1447              LogMessage('warning','Unknown trj version number: %s.' % ICNTRL[0],['gui'])
    1448              return
    1449
    1450          # Diff with doc --> NTRJTI and TRJTIC not in doc
    1451          rec = '!i'
    1452          recSize = struct.calcsize(rec)
    1453          data = struct.unpack(rec, trjFile.read(recSize))
    1454          NTRJTI, = data
    1455
    1456          rec = '!' + '80s'* NTRJTI
    1457          recSize = struct.calcsize(rec)
    1458          trjFile.read(recSize)
    1459
    1460          trjFile.read(calcsize('!8x'))
    1461
    1462          # Record 2
    1463          rec = '!i'
    1464          recSize = struct.calcsize(rec)
    1465          data = struct.unpack(rec, trjFile.read(recSize))
    1466          NEEXTI, = data
    1467
    1468          rec = '!' + '80s'* NEEXTI
    1469          recSize = struct.calcsize(rec)
    1470          trjFile.read(recSize)
    1471
    1472          trjFile.read(calcsize('!8x'))
    1473
    1474          # Record 3
    1475          rec = '!8i8x'
    1476          recSize = struct.calcsize(rec)
    1477          data = struct.unpack(rec, trjFile.read(recSize))
    1478
    1479          PERTYPE, MOLXTL, LCANON, DEFCEL, PRTTHRM, LNOSE, LNPECAN, LTMPDAMP = data
    1480
    1481          LogMessage('info','PERTYPE -- > %s' % PERTYPE,['file'])
    1482          LogMessage('info','MOLXTL -- > %s' % MOLXTL,['file'])
    1483          LogMessage('info','LCANON -- > %s' % LCANON,['file'])
    1484          LogMessage('info','DEFCEL -- > %s' % DEFCEL,['file'])
    1485          LogMessage('info','PRTTHRM -- > %s' % PRTTHRM,['file'])
    1486          LogMessage('info','LNOSE -- > %s' % LNOSE,['file'])
    1487          LogMessage('info','LNPECAN -- > %s' % LNPECAN,['file'])
    1488          LogMessage('info','LTMPDAMP -- > %s' % LTMPDAMP,['file'])
    1489
    1490          # Record 4
    1491          rec = '!i'
    1492          recSize = struct.calcsize(rec)
    1493          data = struct.unpack(rec, trjFile.read(recSize))
    1494          NFLUSD, = data
    1495
    1496          # Diff with doc 8x at the end
    1497          rec = '!%si%si' % (NFLUSD, NFLUSD) + '8s'* NFLUSD + '8x'
    1498          recSize = struct.calcsize(rec)
    1499          data = struct.unpack(rec, trjFile.read(recSize))
    1500
    1501          MVATPF = data[0:NFLUSD]
    1502          NATPFU = data[NFLUSD:2*NFLUSD]
    1503          DECUSD = data[2*NFLUSD:3*NFLUSD]
    1504
    1505          LogMessage('info','MVATPF -- > %s' % str(MVATPF),['file'])
    1506          LogMessage('info','NATPFU -- > %s' % str(NATPFU),['file'])
    1507          LogMessage('info','DECUSD -- > %s' % str(DECUSD),['file'])
    1508
    1509          rec = '!i'
    1510          recSize = struct.calcsize(rec)
    1511          data = struct.unpack(rec, trjFile.read(recSize))
    1512
    1513          TOTMOV = data[0]
    1514
    1515          LogMessage('info','TOTMOV -- > %s' % str(TOTMOV),['file'])
    1516
    1517          rec = '!%si8x' % TOTMOV
    1518          recSize = struct.calcsize(rec)
    1519          data = struct.unpack(rec, trjFile.read(recSize))
    1520
    1521          MVOFST = [d-1 for d in data]
    1522
    1523          LogMessage('info','MVOFST -- > %s' % str(MVOFST),['file'])
    1524
    1525          # Record 4a
    1526          rec = '!i'
    1527          recSize = struct.calcsize(rec)
    1528          data = struct.unpack(rec, trjFile.read(recSize))
    1529
    1530          LEEXTI, = data
    1531
    1532          rec = '!' + '%ss'% LEEXTI
    1533          recSize = struct.calcsize(rec)
    1534          trjFile.read(recSize)
    1535
    1536          trjFile.read(calcsize('!8x'))
    1537
    1538          # Record 4b
    1539          rec = '!i'
    1540          recSize = struct.calcsize(rec)
    1541          data = struct.unpack(rec, trjFile.read(recSize))
    1542
    1543          LPARTI, = data
    1544
    1545          rec = '!' + '%ss'% LPARTI
    1546          recSize = struct.calcsize(rec)
    1547          trjFile.read(recSize)
    1548
    1549          trjFile.read(calcsize('!8x'))
    1550
    1551          while True:
    1552
    1553              try:
    1554
    1555                  # Frame information
    1556                  # Record 1
    1557                  if ICNTRL[0] == 2000:
    1558                      rec = '!%si33%s5i8x' % (prec, prec)
    1559                      velInd = 37
    1560
    1561                  elif ICNTRL[0] == 2010:
    1562                      rec = '!%si57%s6i8x' % (prec, prec)
    1563                      velInd = 61
    1564
    1565                  elif ICNTRL[0] == 3000:
    1566                      rec = '!%si58%s6i8x' % (prec, prec)
    1567                      velInd = 62
    1568
    1569                  recSize = struct.calcsize(rec)
    1570                  data = struct.unpack(rec, trjFile.read(recSize))
    1571
    1572                  CurrentTime     = data[0]
    1573                  itstep          = data[1]
    1574                  Temp            = data[2]
    1575                  AvgTemp         = data[3]
    1576                  TimeStep        = data[4]
    1577                  InitialTemp     = data[5]
    1578                  FinalTemp       = data[6]
    1579                  TotalPE         = data[7]
    1580                  VelocityWritten = data[velInd]
    1581
    1582                  if ICNTRL[0] == 2000:
    1583                      ForcesWritten   = False
    1584
    1585                  elif ICNTRL[0] == 2010:
    1586                      ForcesWritten   = data[62]
    1587                  elif ICNTRL[0] == 3000:
    1588                      ForcesWritten   = data[63]
    1589
    1590                  # Record 2
    1591                  rec = '!12%s8x' % prec
    1592                  recSize = struct.calcsize(rec)
    1593                  data = struct.unpack(rec, trjFile.read(recSize))
    1594
    1595                  Press             = data[0]
    1596                  Volume            = data[1]
    1597                  Junk              = data[2:5]
    1598                  GyrationRadius    = data[5]
    1599                  AvgPress          = data[6]
    1600                  AvgVolume         = data[7]
    1601                  Junk              = data[8:11]
    1602                  AvgGyrationRadius = data[11]
    1603
    1604                  # Record 3
    1605                  if LCANON:
    1606                      rec = '!4%s8x' % prec
    1607                      recSize = struct.calcsize(rec)
    1608                      data = struct.unpack(rec, trjFile.read(recSize))
    1609
    1610                      if LNOSE:
    1611                          snose       = data[0]
    1612                          snoseh      = data[1]
    1613                          dsstot      = data[2]
    1614                          dqcanonNose = data[3]
    1615
    1616                      else:
    1617                          signose   = data[0]
    1618                          zfrict    = data[1]
    1619                          dzprfrict = data[2]
    1620                          dqcanon   = data[3]
    1621
    1622                  # Record 4
    1623                  if DEFCEL > 0:
    1624                      rec = '!22%s8x' % prec
    1625                      recSize = struct.calcsize(rec)
    1626                      data = struct.unpack(rec, trjFile.read(recSize))
    1627
    1628                      a = data[2]*Units.Ang
    1629                      b = data[3]*Units.Ang
    1630                      c = data[4]*Units.Ang
    1631
    1632                      cell = N.zeros((9,), typecode = N.Float)
    1633                      cell[0:3] = a*self.normalizedCell[0:3]
    1634                      cell[3:6] = b*self.normalizedCell[3:6]
    1635                      cell[6:9] = c*self.normalizedCell[6:9]
    1636
    1637                  # Record 5
    1638                  if PERTYPE > 0:
    1639                      rec = '!i14%s8x' % prec
    1640                      recSize = struct.calcsize(rec)
    1641                      data = struct.unpack(rec, trjFile.read(recSize))
    1642
    1643                      NUnitCellAtoms = data[0]
    1644                      Junk           = data[1:15]
    1645
    1646                  # Record 6
    1647                  if LNPECAN:
    1648                      rec = '!3%s8x' % prec
    1649                      recSize = struct.calcsize(rec)
    1650                      trjFile.read(recSize)
    1651
    1652                  # Record 7
    1653                  if LTMPDAMP:
    1654                      rec = '!%s8x' % prec
    1655                      recSize = struct.calcsize(rec)
    1656                      trjFile.read(recSize)
    1657
    1658                  # Record 8
    1659                  rec = '!%s%s' % (TOTMOV, prec)
    1660                  recSize = struct.calcsize(rec)
    1661                  data = struct.unpack(rec, trjFile.read(recSize))
    1662
    1663                  XCOORD = N.array(data[0:TOTMOV])*Units.Ang
    1664                  trjFile.read(calcsize('!8x'))
    1665
    1666                  # Record 9
    1667                  rec = '!%s%s' % (TOTMOV, prec)
    1668                  recSize = struct.calcsize(rec)
    1669                  data = struct.unpack(rec, trjFile.read(recSize))
    1670
    1671                  YCOORD = N.array(data[0:TOTMOV])*Units.Ang
    1672                  trjFile.read(calcsize('!8x'))
    1673
    1674                  # Record 10
    1675                  rec = '!%s%s' % (TOTMOV, prec)
    1676                  recSize = struct.calcsize(rec)
    1677                  data = struct.unpack(rec, trjFile.read(recSize))
    1678
    1679                  ZCOORD = N.array(data[0:TOTMOV])*Units.Ang
    1680                  trjFile.read(calcsize('!8x'))
    1681
    1682                  if VelocityWritten:
    1683                      # Record 11
    1684                      rec = '!%s%s' % (TOTMOV, prec)
    1685                      recSize = struct.calcsize(rec)
    1686                      data = struct.unpack(rec, trjFile.read(recSize))
    1687
    1688                      XVelocity = N.array(data[0:TOTMOV])*Units.Ang
    1689                      trjFile.read(calcsize('!8x'))
    1690
    1691                      # Record 12
    1692                      rec = '!%s%s' % (TOTMOV, prec)
    1693                      recSize = struct.calcsize(rec)
    1694                      data = struct.unpack(rec, trjFile.read(recSize))
    1695
    1696                      YVelocity = N.array(data[0:TOTMOV])*Units.Ang
    1697                      trjFile.read(calcsize('!8x'))
    1698
    1699                      # Record 13
    1700                      rec = '!%s%s' % (TOTMOV, prec)
    1701                      recSize = struct.calcsize(rec)
    1702                      data = struct.unpack(rec, trjFile.read(recSize))
    1703
    1704                      ZVelocity = N.array(data[0:TOTMOV])*Units.Ang
    1705                      trjFile.read(calcsize('!8x'))
    1706
    1707                  if ForcesWritten:
    1708                      # Record 14
    1709                      rec = '!%s%s' % (TOTMOV, prec)
    1710                      recSize = struct.calcsize(rec)
    1711                      trjFile.read(recSize)
    1712
    1713                      trjFile.read(calcsize('!8x'))
    1714
    1715                      # Record 15
    1716                      rec = '!%s%s' % (TOTMOV, prec)
    1717                      recSize = struct.calcsize(rec)
    1718                      trjFile.read(recSize)
    1719
    1720                      trjFile.read(calcsize('!8x'))
    1721
    1722                      # Record 16
    1723                      rec = '!%s%s' % (TOTMOV, prec)
    1724                      recSize = struct.calcsize(rec)
    1725                      trjFile.read(recSize)
    1726
    1727                      trjFile.read(calcsize('!8x'))
    1728
    1729              except:
    1730                  break
    1731
    1732              else:
    1733
    1734                  comp = 0
    1735                  for clust in self.clusters:
    1736
    1737                      for subClust in clust:
    1738                          if subClust['serial_number'] in MVOFST:
    1739                              ind = MVOFST.index(subClust['serial_number'])
    1740                              selCoordinates[comp,0] = XCOORD[ind]
    1741                              selCoordinates[comp,1] = YCOORD[ind]
    1742                              selCoordinates[comp,2] = ZCOORD[ind]
    1743                              selVelocities[comp,0]  = XVelocity[ind]
    1744                              selVelocities[comp,1]  = YVelocity[ind]
    1745                              selVelocities[comp,2]  = ZVelocity[ind]
    1746
    1747                          else:
    1748                              selCoordinates[comp,:] = subClust['coord']
    1749
    1750                          comp += 1
    1751
    1752                  if not self.universe.is_periodic:
    1753                      cell = None
    1754
    1755                  self.universe.setConfiguration(Configuration(self.universe, selCoordinates, cell))
    1756                  self.universe.setVelocities(ParticleVector(self.universe, selVelocities))
    1757
    1758                  snapshot(data = {'time': CurrentTime})
    1759
    1760          trjFile.close()
    1761
    1762          trajectory.close()
    1763
    1764 -    def __convert(self):
    1765          """Performs the actual conversion.
    1766          """
    1767
    1768          self.readXTDFile()
    1769
    1770          if self.module == 'Discover':
    1771              self.readHISFile()
    1772
    1773          elif self.module == 'Forcite':
    1774              self.readTRJFile()
    1775
    1776 -class NAMDConverter:
    1777      """Converts a NAMD Trajectory into a MMTK NetCDFFile.
    1778
    1779      Comments:
    1780          - this code is based on the original converter written by Konrad Hinsen.
    1781      """
    1782
    1783 -    def __init__(self, pdbFile, dcdFile, xstFile, outputFile):
    1784          """
    1785          The constructor. Will do the conversion.
    1786
    1787          @param pdbFile: the NAMD PDB file name of a frame of the trajectory to convert.
    1788          @type pdbFile: string
    1789
    1790          @param dcdFile: the NAMD DCD file name of the trajectory to convert.
    1791          @type dcdFile: string
    1792
    1793          @param xstFile: the NAMD XSTfile name of the trajectory to convert.
    1794          @type xstFile: string
    1795
    1796          @param outputFile: the name of MMTK NetCDF trajectory output file.
    1797          @type outputFile: string
    1798          """
    1799
    1800          self.pdbFile = pdbFile
    1801          self.dcdFile = dcdFile
    1802          self.xstFile = xstFile
    1803          self.outputFile = outputFile
    1804
    1805          # Do the conversion.
    1806          self.__convert()
    1807
    1808 -    def __convert(self):
    1809          """Performs the actual conversion."""
    1810
    1811          # First process the XST file if one was loaded.
    1812          if self.xstFile is not None:
    1813              xst = open(self.xstFile, 'r')
    1814              data = xst.readlines()[3:]
    1815              xst.close()
    1816              xstData = []
    1817              for line in data:
    1818                  d = N.reshape(N.array([float(v) for v in line.split()[1:10]]),(3,3))
    1819                  xstData.append([Vector(v)*Units.Ang for v in d])
    1820
    1821              universe = ParallelepipedicPeriodicUniverse(xstData[0])
    1822
    1823          else:
    1824              universe = InfiniteUniverse()
    1825
    1826          # Open the DCD trajectory file
    1827          dcd = DCDFile(self.dcdFile)
    1828          step = dcd.istart
    1829          step_increment = dcd.nsavc
    1830          dt = dcd.delta
    1831          unit_cell, x, y, z = dcd.readStep()
    1832
    1833          # Create all objects from the PDB file. The PDB file must match the
    1834          # the DCD file (same atom order).
    1835          conf = PDBConfiguration(self.pdbFile)
    1836          universe.addObject(conf.createAll())
    1837
    1838          # Create the output trajectory
    1839          trajectory = Trajectory(universe, self.outputFile, mode = 'w', comment = dcd.title)
    1840          snapshot = SnapshotGenerator(universe, actions=[TrajectoryOutput(trajectory, ["all"], 0, None, 1)])
    1841
    1842          # Iterate over the DCD file and write to the output trajectory
    1843          conf_array = universe.configuration().array
    1844          comp = 0
    1845          while True:
    1846              conf_array[:, 0] = x
    1847              conf_array[:, 1] = y
    1848              conf_array[:, 2] = z
    1849
    1850              if universe.is_periodic:
    1851                  universe.setShape(xstData[comp])
    1852
    1853              step_data = {'time': step*dt, 'step': step}
    1854              snapshot(data = step_data)
    1855              step += step_increment
    1856
    1857              try:
    1858                  unit_cell, x, y, z = dcd.readStep()
    1859                  comp += 1
    1860
    1861              except EndOfFile:
    1862                  break
    1863
    1864          # Close the output trajectory
    1865          trajectory.close()
    1866
    1867 -class VASPConverter:
    1868      """Converts a VASP Trajectory into a MMTK NetCDFFile."""
    1869
    1870 -    def __init__(self, contcarFile, xdatcarFile, outputFile, atomContents):
    1871          """
    1872          The constructor. Will do the conversion.
    1873
    1874          @param contcarFile: the VASP CONTCAR or POSCAR file name of the trajectory to convert.
    1875          @type contcarFile: string
    1876
    1877          @param xdatcarFile: the VASP XDATCAR file name of the trajectory to convert.
    1878          @type xdatcarFile: string
    1879
    1880          @param outputFile: the name of MMTK NetCDF trajectory output file.
    1881          @type outputFile: string
    1882
    1883          @param atomContents: List of the element names (string) in the order they appear in the
    1884              trajectory.
    1885          @type atomContents: list
    1886          """
    1887
    1888          self.contcarFile = contcarFile
    1889          self.xdatcarFile = xdatcarFile
    1890          self.outputFile = outputFile
    1891          self.atomContents = atomContents
    1892
    1893          # Do the conversion.
    1894          self.__convert()
    1895
    1896 -    def __convert(self):
    1897          """Performs the actual conversion."""
    1898
    1899          # The VASP/CONTCAR is opened for reading.
    1900          contcarFile = open(self.contcarFile, 'r')
    1901
    1902          # The first line is read. It may contains the name of the system.
    1903          firstLine = contcarFile.readline().split()
    1904          # Guess the name of the system from the first line of the CONTCAR.
    1905          try:
    1906              clustName = firstLine[firstLine.index('System=')+1]
    1907          except:
    1908              clustName = 'Cluster'
    1909
    1910          # Pick up the scale factor that is the first element of the second line of CONTCAR file.
    1911          scale_factor = float(contcarFile.readline().split()[0])
    1912
    1913          # The direct basis vectors.
    1914          dirA = Vector([float(v) for v in contcarFile.readline().split()])*Units.Ang*scale_factor
    1915          dirB = Vector([float(v) for v in contcarFile.readline().split()])*Units.Ang*scale_factor
    1916          dirC = Vector([float(v) for v in contcarFile.readline().split()])*Units.Ang*scale_factor
    1917
    1918          # In VASP, there are always some PBC.
    1919          universe = ParallelepipedicPeriodicUniverse((dirA, dirB, dirC), None)
    1920
    1921          # The atomic composition of the molecule. A list of integer giving how many number of atom of each type is present
    1922          # in the molecule
    1923          numberOfAtomsPerType = [int(v) for v in contcarFile.readline().split()]
    1924
    1925          # Close |CONTCAR| VASP file right now. No need for it anymore.
    1926          contcarFile.close()
    1927
    1928          mmtkAtomList = []
    1929          # The atoms of the system are added to the universe.
    1930          for i in range(len(self.atomContents)):
    1931              for l in range(numberOfAtomsPerType[i]):
    1932                  mmtkAtomList.append(Atom(self.atomContents[i], name = self.atomContents[i] + str(l+1)))
    1933
    1934          universe.addObject(AtomCluster(mmtkAtomList, name = clustName))
    1935
    1936          # Get data from xdatcar file.
    1937          xdatcarFile = open(self.xdatcarFile, 'r')
    1938          nAtoms = int(xdatcarFile.readline().split()[0])
    1939
    1940          if nAtoms != universe.numberOfAtoms():
    1941              raise Error('Wrong number of atoms.')
    1942
    1943          # The last entry of the second line of the XDATCAR file gives the timestep in seconds. Its is converted in ps.
    1944          timestep = float(xdatcarFile.readline().split()[-1])*Units.s
    1945          xdatcarFile.readline()
    1946          xdatcarFile.readline()
    1947          xdatcarFile.readline()
    1948
    1949          # The output trajectory NetCDF file is opened for writing.
    1950          trajectory = Trajectory(universe, self.outputFile, mode = 'w', comment = 'From VASP file')
    1951          snapshot = SnapshotGenerator(universe, actions = [TrajectoryOutput(trajectory,["all"], 0, None, 1)])
    1952
    1953          step = 0
    1954          while True:
    1955              t = step*timestep
    1956              # This is a blank line separating each frame.
    1957              junk = xdatcarFile.readline()
    1958
    1959              if not junk:
    1960                  break
    1961
    1962              try:
    1963                  # Loop over the coordinates for the current frame.
    1964                  for at in range(len(universe.atomList())):
    1965                      line = xdatcarFile.readline()
    1966                      # Conversion of the coordinates into a Scientific Vector object
    1967                      coord = Vector([float(v) for v in line.split()[0:3]])
    1968                      # In VASP the coordinates are given in box coordinates. They must be converted in real coordinates.
    1969                      coord = universe.boxToRealCoordinates(coord)
    1970                      # The universe is updated with the current position of atom |at|.
    1971                      universe.atomList()[at].setPosition(coord)
    1972                  # The whole frame plus the |time| extra variable are flushed into the output file.
    1973                  snapshot(data = {'time': t})
    1974                  step += 1
    1975
    1976              except:
    1977                  break
    1978
    1979          # The input and output file are closed.
    1980          xdatcarFile.close()
    1981          trajectory.close()
    1982