#!/usr/bin/python
# Find the top of the valence band and the bottom of the conduction band from EIGENVAL  
# Usage: chmod +x gbandedges
#        gbandedges [path] [ef]
#        ef: trial Fermi energy
#        ef read from OUTCAR if not given
# Last modified on Jun.21,2014

from sys import argv,exit
from os import path

ief = 0
if len(argv)==1:
    dir = '.'
elif len(argv)==2:
    try:
        ef = float(argv[1])
        dir = '.'
        ief = 1
    except: 
        dir = argv[1].rstrip('/')
elif len(argv)==3:
    dir = argv[1].rstrip('/')
    ef = float(argv[2])
    ief = 1
else:
    print "Usage: gbandedges [path] [Ef]"
    print "       ef: trial Fermi energy"
    print "       ef read from OUTCAR if not given"
    exit(0)

fn_outcar = dir+'/OUTCAR'
if not ief:
    if path.exists(fn_outcar):
        for line in open(fn_outcar, 'r'):
            if line.find('E-fermi') > -1:
                ef=float(line.split()[2])
    else:
        print 'OUTCAR does not exist!'    
        exit(0)

fn_eig = dir+'/EIGENVAL'
if not path.exists(fn_eig):
    print 'EIGENVAL does not exist!'    
    exit(0)
    
f_eig=open(fn_eig)
lines=f_eig.readlines()
file.close
nspin=int(lines[0].split()[3])
nkpoint=int(lines[5].split()[1])
nband=int(lines[5].split()[2])

eup=ef+99.0
edown=ef-99.0
k=[]
for i in range(nkpoint):
    k.append(lines[7+i*(nband+2)].split()[0:3])
    for j in range(nband):
        for s in range(nspin):
            eval=float(lines[8+i*(nband+2)+j].split()[s+1])
            d=eval-ef 
            if (d>0):
                if (eval<eup):
                    eup=eval
                    kup=i
                    bup=j
                    sup=s
            if (d<0):
                if (eval>edown):
                    edown=eval
                    kdown=i
                    bdown=j
                    sdown=s

print 'Efermi  : %.3f'%ef
print 'N band  :',nband
print 'N kpoint:',nkpoint
print 'N spin  :',nspin
print 'Band gap: %.3f'%(eup-edown)
print 'Bottom of conduction band:'
print '  Energy : %.3f'%eup
print '  K point:',kup+1,k[kup]
print '  Band   :',bup+1
print '  Spin   :',sup+1
print 'Top of valence band:'
print '  Energy : %.3f'%edown
print '  K point:',kdown+1,k[kdown]
print '  Band   :',bdown+1
print '  Spin   :',sdown+1