qbox_dos.py 2.71 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 [-last] 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
# With the -last option, only the last <eigenset> is used to compute the DOS
Francois Gygi committed
9 10 11 12 13

import xml.sax
import sys
import math

Francois Gygi committed
14 15
if (len(sys.argv) != 5) and (len(sys.argv) != 6) :
  print "use: ",sys.argv[0]," [-last] emin emax width file.r"
Francois Gygi committed
16 17
  sys.exit()

Francois Gygi committed
18 19 20 21 22 23 24 25 26 27 28 29 30
iarg = 1
lastonly = False
if (sys.argv[iarg] == "-last") :
  lastonly = True
  iarg += 1

emin = float(sys.argv[iarg])
iarg += 1
emax = float(sys.argv[iarg])
iarg += 1
width = float(sys.argv[iarg])
iarg += 1
infile = sys.argv[iarg]
Francois Gygi committed
31 32 33 34 35 36 37 38 39 40 41 42

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):
43
    self.nspin = 1
Francois Gygi committed
44
    self.readData = 0
45 46
    self.dos_up = [0] * ndos
    self.dos_dn = [0] * ndos
Francois Gygi committed
47 48

  def startElement(self, name, attributes):
Francois Gygi committed
49 50 51
    if (name == "eigenset") and (lastonly):
      self.dos_up = [0] * ndos
      self.dos_dn = [0] * ndos
Francois Gygi committed
52 53
    if name == "eigenvalues":
      self.n = attributes["n"]
54
      self.spin = int(attributes["spin"])
Francois Gygi committed
55
      self.kpoint = attributes["kpoint"]
56
      self.weight = float(attributes["weight"])
Francois Gygi committed
57 58
      self.readData = 1
      self.buffer = ""
59 60
      if self.spin == 1:
        self.nspin = 2
Francois Gygi committed
61 62 63 64 65 66 67 68

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

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

71
  def accumulate_dos(self):
Francois Gygi committed
72
    self.e = self.buffer.split()
73 74 75 76 77 78 79 80 81 82 83 84 85
    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
86 87
    for j in range(ndos):
      ej = emin + j * de
88 89 90 91 92 93 94 95
      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
96 97 98 99 100

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