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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
|
import functools
import warnings
import numpy as np
from ase.utils import IOContext
def get_band_gap(calc, direct=False, spin=None, output='-'):
warnings.warn('Please use ase.dft.bandgap.bandgap() instead!')
gap, (s1, k1, n1), (s2, k2, n2) = bandgap(calc, direct, spin, output)
ns = calc.get_number_of_spins()
if ns == 2 and spin is None:
return gap, (s1, k1), (s2, k2)
return gap, k1, k2
def bandgap(calc=None, direct=False, spin=None, output='-',
eigenvalues=None, efermi=None, kpts=None):
"""Calculates the band-gap.
Parameters:
calc: Calculator object
Electronic structure calculator object.
direct: bool
Calculate direct band-gap.
spin: int or None
For spin-polarized systems, you can use spin=0 or spin=1 to look only
at a single spin-channel.
output: file descriptor
Use output=None for no text output or '-' for stdout (default).
eigenvalues: ndarray of shape (nspin, nkpt, nband) or (nkpt, nband)
Eigenvalues.
efermi: float
Fermi level (defaults to 0.0).
kpts: ndarray of shape (nkpt, 3)
For pretty text output only.
Returns a (gap, p1, p2) tuple where p1 and p2 are tuples of indices of the
valence and conduction points (s, k, n).
Example:
>>> gap, p1, p2 = bandgap(silicon.calc)
Gap: 1.2 eV
Transition (v -> c):
[0.000, 0.000, 0.000] -> [0.500, 0.500, 0.000]
>>> print(gap, p1, p2)
1.2 (0, 0, 3), (0, 5, 4)
>>> gap, p1, p2 = bandgap(silicon.calc, direct=True)
Direct gap: 3.4 eV
Transition at: [0.000, 0.000, 0.000]
>>> print(gap, p1, p2)
3.4 (0, 0, 3), (0, 0, 4)
"""
if calc:
kpts = calc.get_ibz_k_points()
nk = len(kpts)
ns = calc.get_number_of_spins()
eigenvalues = np.array([[calc.get_eigenvalues(kpt=k, spin=s)
for k in range(nk)]
for s in range(ns)])
if efermi is None:
efermi = calc.get_fermi_level()
efermi = efermi or 0.0
e_skn = eigenvalues - efermi
if eigenvalues.ndim == 2:
e_skn = e_skn[np.newaxis] # spinors
if not np.isfinite(e_skn).all():
raise ValueError('Bad eigenvalues!')
gap, (s1, k1, n1), (s2, k2, n2) = _bandgap(e_skn, spin, direct)
with IOContext() as iocontext:
fd = iocontext.openfile(output)
p = functools.partial(print, file=fd)
def skn(s, k, n):
"""Convert k or (s, k) to string."""
if kpts is None:
return '(s={}, k={}, n={})'.format(s, k, n)
return '(s={}, k={}, n={}, [{:.2f}, {:.2f}, {:.2f}])'.format(
s, k, n, *kpts[k])
if spin is not None:
p('spin={}: '.format(spin), end='')
if gap == 0.0:
p('No gap')
elif direct:
p('Direct gap: {:.3f} eV'.format(gap))
if s1 == s2:
p('Transition at:', skn(s1, k1, n1))
else:
p('Transition at:', skn('{}->{}'.format(s1, s2), k1, n1))
else:
p('Gap: {:.3f} eV'.format(gap))
p('Transition (v -> c):')
p(' ', skn(s1, k1, n1), '->', skn(s2, k2, n2))
if eigenvalues.ndim != 3:
p1 = (k1, n1)
p2 = (k2, n2)
else:
p1 = (s1, k1, n1)
p2 = (s2, k2, n2)
return gap, p1, p2
def _bandgap(e_skn, spin, direct):
"""Helper function."""
ns, nk, nb = e_skn.shape
s1 = s2 = k1 = k2 = n1 = n2 = None
N_sk = (e_skn < 0.0).sum(2) # number of occupied bands
# Check for bands crossing the fermi-level
if ns == 1:
if N_sk[0].ptp() > 0:
return 0.0, (None, None, None), (None, None, None)
elif spin is None:
if (N_sk.ptp(axis=1) > 0).any():
return 0.0, (None, None, None), (None, None, None)
elif N_sk[spin].ptp() > 0:
return 0.0, (None, None, None), (None, None, None)
if (N_sk == 0).any() or (N_sk == nb).any():
raise ValueError('Too few bands!')
e_skn = np.array([[e_skn[s, k, N_sk[s, k] - 1:N_sk[s, k] + 1]
for k in range(nk)]
for s in range(ns)])
ev_sk = e_skn[:, :, 0] # valence band
ec_sk = e_skn[:, :, 1] # conduction band
if ns == 1:
s1 = 0
s2 = 0
gap, k1, k2 = find_gap(ev_sk[0], ec_sk[0], direct)
n1 = N_sk[0, 0] - 1
n2 = n1 + 1
return gap, (0, k1, n1), (0, k2, n2)
if spin is None:
gap, k1, k2 = find_gap(ev_sk.ravel(), ec_sk.ravel(), direct)
if direct:
# Check also spin flips:
for s in [0, 1]:
g, k, _ = find_gap(ev_sk[s], ec_sk[1 - s], direct)
if g < gap:
gap = g
k1 = k + nk * s
k2 = k + nk * (1 - s)
if gap > 0.0:
s1, k1 = divmod(k1, nk)
s2, k2 = divmod(k2, nk)
n1 = N_sk[s1, k1] - 1
n2 = N_sk[s2, k2]
return gap, (s1, k1, n1), (s2, k2, n2)
return 0.0, (None, None, None), (None, None, None)
gap, k1, k2 = find_gap(ev_sk[spin], ec_sk[spin], direct)
s1 = spin
s2 = spin
n1 = N_sk[s1, k1] - 1
n2 = n1 + 1
return gap, (s1, k1, n1), (s2, k2, n2)
def find_gap(ev_k, ec_k, direct):
"""Helper function."""
if direct:
gap_k = ec_k - ev_k
k = gap_k.argmin()
return gap_k[k], k, k
kv = ev_k.argmax()
kc = ec_k.argmin()
return ec_k[kc] - ev_k[kv], kv, kc
|
