qbox_dos.py 2.36 KB
Newer Older
Francois Gygi committed
1
#!/usr/bin/python
2 3
# qbox_dos.py: extract electronic DOS from Qbox output
# generate DOS plot in gnuplot format
Francois Gygi committed
4
# use: qbox_dos.py emin emax width file.r
5 6
# emin, emax: bounds of plot in [eV]
# width: gaussian broadening in [eV]
7
# the DOS is accumulated separately for each spin
Francois Gygi committed
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

import xml.sax
import sys
import math

if len(sys.argv) != 5:
  print "use: ",sys.argv[0]," emin emax width file.r"
  sys.exit()

emin = float(sys.argv[1])
emax = float(sys.argv[2])
width = float(sys.argv[3])
infile = sys.argv[4]

ndos = 501
de = (emax - emin)/(ndos-1)

# normalized gaussian distribution in one dimension
# f(x) = 1/(sqrt(pi)*width) * exp(-(x/width)^2 )
def gauss(x, width):
  return (1.0/(math.sqrt(math.pi)*width)) * math.exp(-(x/width)**2)
  
# Qbox output handler to extract and process data
class QboxOutputHandler(xml.sax.handler.ContentHandler):
  def __init__(self):
33
    self.nspin = 1
Francois Gygi committed
34
    self.readData = 0
35 36
    self.dos_up = [0] * ndos
    self.dos_dn = [0] * ndos
Francois Gygi committed
37 38 39 40

  def startElement(self, name, attributes):
    if name == "eigenvalues":
      self.n = attributes["n"]
41
      self.spin = int(attributes["spin"])
Francois Gygi committed
42
      self.kpoint = attributes["kpoint"]
43
      self.weight = float(attributes["weight"])
Francois Gygi committed
44 45
      self.readData = 1
      self.buffer = ""
46 47
      if self.spin == 1:
        self.nspin = 2
Francois Gygi committed
48 49 50 51 52 53 54 55

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

  def endElement(self, name):
    if name == "eigenvalues":
      self.readData = 0
56
      self.accumulate_dos()
Francois Gygi committed
57

58
  def accumulate_dos(self):
Francois Gygi committed
59
    self.e = self.buffer.split()
60 61 62 63 64 65 66 67 68 69 70 71 72
    if self.spin == 0:
      for i in range(len(self.e)):
        for j in range(ndos):
          ej = emin + j * de
          self.dos_up[j] += gauss(float(self.e[i])-ej, width ) * self.weight
    if self.spin == 1:
      for i in range(len(self.e)):
        for j in range(ndos):
          ej = emin + j * de
          self.dos_dn[j] += gauss(float(self.e[i])-ej, width ) * self.weight

  def print_dos(self):
    print "# ",infile," spin=0 width=",width
Francois Gygi committed
73 74
    for j in range(ndos):
      ej = emin + j * de
75 76 77 78 79 80 81 82
      print ej, self.dos_up[j]
    if self.nspin == 2:
      print
      print
      print "# ",infile," spin=1 width=",width
      for j in range(ndos):
        ej = emin + j * de
        print ej, self.dos_dn[j]
Francois Gygi committed
83 84 85 86 87

parser = xml.sax.make_parser()
handler = QboxOutputHandler()
parser.setContentHandler(handler)
parser.parse(infile)
88
handler.print_dos()