1 #!/usr/bin/env python
   2 """Program to plot the density of states for a Dacapo calculation
   3 
   4 Make a gnuplot:
   5 
   6 > dos.py filename.nc
   7 
   8 Dump the data to a file:
   9 
  10 > dos.py filename.nc dos.dat
  11 """
  12 
  13 import sys
  14 from Dacapo import Dacapo
  15 from ASE.Utilities.DOS import DOS
  16 import Gnuplot as gp
  17 
  18 # Get the filename of the Dacapo output file:
  19 filename = sys.argv[1]
  20 
  21 # Read in the atom and restart the calculator:
  22 atoms = Dacapo.ReadAtoms(filename)
  23 
  24 # Create a DOS object and extract data:
  25 dos = DOS(atoms, width=0.2)
  26 e = dos.GetEnergies()
  27 dos0 = dos.GetDOS(0)
  28 if dos.nspins == 2:
  29     # Spin polarized calculation:
  30     dos1 = dos.GetDOS(1)
  31 
  32 if len(sys.argv) > 2:
  33     # Dump the numbers to a file:
  34     out = file(sys.argv[2], 'w')
  35     if dos.nspins == 1:
  36         for x, y in zip(e, dos0):
  37             print >> out, x, y
  38     else:
  39         for x, y1, y2 in zip(e, dos0, dos1):
  40             print >> out, x, y1, y2
  41 else:
  42     # Make a gnuplot of the data:
  43     plot = gp.Gnuplot(persist=True)
  44     if dos.nspins == 1:
  45         plot.plot(gp.Data(e, dos0, with='lines'))
  46     else:
  47         plot.plot(gp.Data(e, dos0, with='lines'),
  48                   gp.Data(e, dos1, with='lines'))