qbox_msd.py 1.75 KB
Newer Older
1
#!/usr/bin/python
2 3 4
# qbox_msd.py: compute mean-square displacement in an MD simulation
# generate plot of <r^2>(t) in gnuplot format
# use: qbox_msd.py species file.r [file.r ...]
5 6 7 8 9

import xml.sax
import sys
import math

10 11
if len(sys.argv) < 3:
  print "use: ",sys.argv[0]," species file [file ...]"
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
  sys.exit()

species = sys.argv[1]

# Qbox output handler to extract and process data
class QboxOutputHandler(xml.sax.handler.ContentHandler):
  def __init__(self):
    self.step = 0
    self.readData = 0
    self.readAtom = 0
    self.tau0 = []

  def startElement(self, name, attributes):
    if name == "atomset":
      self.tau=[]
    elif name == "atom":
      self.readAtom = (attributes["species"] == species)
    elif (name == "position") & self.readAtom:
        self.readData = 1
        self.buffer = ""

  def characters(self, data):
    if self.readData:
      self.buffer += data

  def endElement(self, name):
    if (name == "position") & self.readAtom:
      self.readData = 0
      pos = self.buffer.split()
      x = float(pos[0])
      y = float(pos[1])
      z = float(pos[2])
      self.tau.append([x,y,z])
    elif name == "atomset":
      if self.step == 0:
        # copy initial positions to tau0
        self.tau0 = self.tau
      # compute square displacement
      disp2sum = 0.0
      for i in range(len(self.tau)):
        dx = self.tau[i][0]-self.tau0[i][0]
        dy = self.tau[i][1]-self.tau0[i][1]
        dz = self.tau[i][2]-self.tau0[i][2]
        disp2sum += dx*dx + dy*dy + dz*dz
56
      print '%12.6f'%(disp2sum/len(self.tau))
57 58
      self.step += 1

59
print "# ",species
60 61 62
parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
63 64 65 66
for i in range(len(sys.argv)-2):
  infile = sys.argv[i+2]
  parser.parse(infile)

67 68
print
print