forked from yyyu200/QEbandplot
-
Notifications
You must be signed in to change notification settings - Fork 0
/
bandplot.py
114 lines (86 loc) · 2.95 KB
/
bandplot.py
1
2
3
4
5
6
7
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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
# -*- coding: utf-8 -*-
"""
Author: yyyu200@163.com
Modified by: hongyi.zhao@gmail.com
"""
import os
import re
import matplotlib.pyplot as plt
import numpy as np
# https://www.quantum-espresso.org/Doc/INPUT_BANDS.html#idm19
# filband CHARACTER
#Default: 'bands.out'
# file name for band output (to be read by "plotband.x")
def parse_filband(feig, npl=10):
# feig : filband in bands.x input file
# npl : number per line, 10 for bands.x, 6 for phonon
f = open(feig, 'r')
lines = f.readlines()
header = lines[0].strip()
line = header.strip('\n')
shape = re.split('[,=/]', line)
nbnd = int(shape[1])
nks = int(shape[3])
eig = np.zeros((nks, nbnd+1), dtype=np.float32)
dividend = nbnd
divisor = npl
div = nbnd // npl + 1 if nbnd % npl == 0 else nbnd // npl + 2
kinfo = []
for index, value in enumerate(lines[1:]):
value = value.strip(' \n')
quotient = index // div
remainder = index % div
if remainder == 0:
kinfo.append(value)
else:
value = re.split('[ ]+', value)
a = (remainder - 1) * npl
b = a + len(value)
eig[quotient][a:b] = value
f.close()
return eig, nbnd, nks, kinfo
do_find_gap = True
e_ref = 0.0 # set to fermi-level in scf output for metal, only applicable for do_find_gap=False
if do_find_gap:
nvband = 26 # valence band number, only applicable for do_find_gap=True
else:
nvband = 0
ymin = -10 # y range in plot
ymax = 10
lw = 1.2 # line width
p1 = plt.subplot(1, 1, 1)
F = plt.gcf()
# F.set_size_inches([5,5])
eig, nbnd, nks, kinfo = parse_filband('bands.out')
if nbnd <= nvband:
print("warning: nvband should be less than the calculated band.")
plt.xlim([0, nks-1]) # k-points
plt.ylim([ymin, ymax])
#plt.xlabel(r'$k (\AA^{-1})$',fontsize=16)
plt.ylabel(r' E (eV) ', fontsize=16)
# for insulators only, nvband can be found by gappw.sh(https://github.com/yyyu200/gappw)
if do_find_gap and nbnd > nvband:
eig_vbm = max(eig[:, nvband-1])
eig_cbm = min(eig[:, nvband])
Gap = eig_cbm-eig_vbm
plt.title("Band gap= %.4f eV" % (Gap))
e_ref = eig_vbm
for i in range(nbnd):
line1 = plt.plot(eig[:, i]-e_ref, color='r', linewidth=lw)
# positions of vertical lines, or specified by [0, 20, 40, ...]
vlines = np.arange(0, nks, 20)
for vline in vlines:
plt.axvline(x=vline, ymin=ymin, ymax=ymax, linewidth=lw, color='black')
xlabeltext = [r'${\Gamma}$', 'X', 'M', r'${\Gamma}$',
'Z', 'R', 'A', 'Z', 'X', 'R', 'M', 'A']
if len(xlabeltext) < len(vlines):
for i in range(len(vlines)-len(xlabeltext)):
xlabeltext.append('X')
elif len(xlabeltext) > len(vlines):
xlabeltext = xlabeltext[0:len(vlines)]
assert len(xlabeltext)==len(vlines)
plt.xticks( vlines, xlabeltext)
plt.text(4, 8, 'SnO$_{2}$', fontsize=12, color='black', bbox=dict(facecolor='white',alpha=0.99,edgecolor='black') )
plt.tight_layout()
plt.savefig('pwband.png',dpi=500)
plt.show()