diff --git a/util/qbox_dos/qbox_dos.py b/util/qbox_dos.py
similarity index 96%
rename from util/qbox_dos/qbox_dos.py
rename to util/qbox_dos.py
index d0d4059..766373b 100755
--- a/util/qbox_dos/qbox_dos.py
+++ b/util/qbox_dos.py
@@ -1,7 +1,9 @@
 #!/usr/bin/python
-# qbox_dos.py: extract dos from Qbox output
-# generate dos plot in gnuplot format
+# qbox_dos.py: extract electronic DOS from Qbox output
+# generate DOS plot in gnuplot format
 # use: qbox_dos.py emin emax width file.r
+# emin, emax: bounds of plot in [eV]
+# width: gaussian broadening in [eV]
 # the DOS is accumulated separately for each spin
 
 import xml.sax
diff --git a/util/qbox_msd.py b/util/qbox_msd.py
new file mode 100755
index 0000000..3f7f685
--- /dev/null
+++ b/util/qbox_msd.py
@@ -0,0 +1,68 @@
+#!/usr/bin/python
+# 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 ...]
+
+import xml.sax
+import sys
+import math
+
+if len(sys.argv) < 3:
+  print "use: ",sys.argv[0]," species file [file ...]"
+  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
+      print self.step,'%12.6f'%(disp2sum/len(self.tau))
+      self.step += 1
+
+print "# ",species
+parser = xml.sax.make_parser()
+handler = QboxOutputHandler()
+parser.setContentHandler(handler)
+for i in range(len(sys.argv)-2):
+  infile = sys.argv[i+2]
+  parser.parse(infile)
+
+print
+print
diff --git a/util/qbox_repair_h2o.py b/util/qbox_repair_h2o.py
new file mode 100755
index 0000000..dfdcb0a
--- /dev/null
+++ b/util/qbox_repair_h2o.py
@@ -0,0 +1,58 @@
+#!/usr/bin/python
+# qbox_repair_h2o.py: repair broken h2o molecules in a Qbox sys file
+# move hydrogen atoms across periodic boundaries to repair molecules
+# use: qbox_repair_h2o.py file.sys
+
+import sys
+import math
+
+olist = []
+hlist = []
+
+def distance(a,b,sx,sy,sz):
+  return math.sqrt((a[3]-b[3]-sx)**2+(a[4]-b[4]-sy)**2+(a[5]-b[5]-sz)**2)
+
+f = open(sys.argv[1])
+for line in f:
+  l = line.split()
+  if ( l[0] == "set" ) & ( l[1] == "cell" ):
+    a_cell = float(l[2])
+    b_cell = float(l[6])
+    c_cell = float(l[10])
+    print line
+  elif ( l[0] == "species" ):
+    print line
+  elif ( l[0] == "atom" ) & ( l[2] == "oxygen" ):
+    olist.append([l[0],l[1],l[2],float(l[3]),float(l[4]),float(l[5])])
+  elif ( l[0] == "atom" ) & ( l[2] == "hydrogen" ):
+    hlist.append([l[0],l[1],l[2],float(l[3]),float(l[4]),float(l[5])])
+
+for h in hlist:
+  # find nearest oxygen atom in olist
+  mindist = 1.e10
+  for o in olist:
+    # compute minimal distance o-h
+    for sx in [-a_cell,0,a_cell]:
+      for sy in [-b_cell,0,b_cell]:
+        for sz in [-c_cell,0,c_cell]:
+          d = distance(o,h,sx,sy,sz)
+          # print "dist(",sx,sy,sz,") = ",d
+          if d < mindist:
+            mindist = d
+            sx_min = sx
+            sy_min = sy
+            sz_min = sz
+            nearest_o = o;
+  # print "min shift is: ",sx_min,sy_min,sz_min
+  if ( sx_min != 0 ) | ( sy_min != 0 ) | ( sz_min != 0 ):
+    print "# current ",h[1]," at ", h[3],h[4],h[5]
+    print "# nearest O is at ", nearest_o[3],nearest_o[4],nearest_o[5]
+    print "# move ",h[1]," by ", sx_min, sy_min, sz_min
+    h[3] += sx_min    
+    h[4] += sy_min    
+    h[5] += sz_min    
+
+for o in olist:
+  print o[0],o[1],o[2],'%10.5f'%o[3],'%10.5f'%o[4],'%10.5f'%o[5]
+for h in hlist:
+  print h[0],h[1],h[2],'%10.5f'%h[3],'%10.5f'%h[4],'%10.5f'%h[5]