r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1import functools

2import warnings

4import numpy as np

6from ase.utils import IOContext

9def get_band_gap(calc, direct=False, spin=None, output='-'):

11 gap, (s1, k1, n1), (s2, k2, n2) = bandgap(calc, direct, spin, output)

12 ns = calc.get_number_of_spins()

13 if ns == 2 and spin is None:

14 return gap, (s1, k1), (s2, k2)

15 return gap, k1, k2

18def bandgap(calc=None, direct=False, spin=None, output='-',

19 eigenvalues=None, efermi=None, kpts=None):

20 """Calculates the band-gap.

22 Parameters:

24 calc: Calculator object

25 Electronic structure calculator object.

26 direct: bool

27 Calculate direct band-gap.

28 spin: int or None

29 For spin-polarized systems, you can use spin=0 or spin=1 to look only

30 at a single spin-channel.

31 output: file descriptor

32 Use output=None for no text output or '-' for stdout (default).

33 eigenvalues: ndarray of shape (nspin, nkpt, nband) or (nkpt, nband)

34 Eigenvalues.

35 efermi: float

36 Fermi level (defaults to 0.0).

37 kpts: ndarray of shape (nkpt, 3)

38 For pretty text output only.

40 Returns a (gap, p1, p2) tuple where p1 and p2 are tuples of indices of the

41 valence and conduction points (s, k, n).

43 Example:

45 >>> gap, p1, p2 = bandgap(silicon.calc)

46 Gap: 1.2 eV

47 Transition (v -> c):

48 [0.000, 0.000, 0.000] -> [0.500, 0.500, 0.000]

49 >>> print(gap, p1, p2)

50 1.2 (0, 0, 3), (0, 5, 4)

51 >>> gap, p1, p2 = bandgap(silicon.calc, direct=True)

52 Direct gap: 3.4 eV

53 Transition at: [0.000, 0.000, 0.000]

54 >>> print(gap, p1, p2)

55 3.4 (0, 0, 3), (0, 0, 4)

56 """

58 if calc:

59 kpts = calc.get_ibz_k_points()

60 nk = len(kpts)

61 ns = calc.get_number_of_spins()

62 eigenvalues = np.array([[calc.get_eigenvalues(kpt=k, spin=s)

63 for k in range(nk)]

64 for s in range(ns)])

65 if efermi is None:

66 efermi = calc.get_fermi_level()

68 efermi = efermi or 0.0

70 e_skn = eigenvalues - efermi

71 if eigenvalues.ndim == 2:

72 e_skn = e_skn[np.newaxis] # spinors

74 if not np.isfinite(e_skn).all():

77 gap, (s1, k1, n1), (s2, k2, n2) = _bandgap(e_skn, spin, direct)

79 with IOContext() as iocontext:

80 fd = iocontext.openfile(output)

81 p = functools.partial(print, file=fd)

83 def skn(s, k, n):

84 """Convert k or (s, k) to string."""

85 if kpts is None:

86 return '(s={}, k={}, n={})'.format(s, k, n)

87 return '(s={}, k={}, n={}, [{:.2f}, {:.2f}, {:.2f}])'.format(

88 s, k, n, *kpts[k])

90 if spin is not None:

91 p('spin={}: '.format(spin), end='')

92 if gap == 0.0:

93 p('No gap')

94 elif direct:

95 p('Direct gap: {:.3f} eV'.format(gap))

96 if s1 == s2:

97 p('Transition at:', skn(s1, k1, n1))

98 else:

99 p('Transition at:', skn('{}->{}'.format(s1, s2), k1, n1))

100 else:

101 p('Gap: {:.3f} eV'.format(gap))

102 p('Transition (v -> c):')

103 p(' ', skn(s1, k1, n1), '->', skn(s2, k2, n2))

105 if eigenvalues.ndim != 3:

106 p1 = (k1, n1)

107 p2 = (k2, n2)

108 else:

109 p1 = (s1, k1, n1)

110 p2 = (s2, k2, n2)

112 return gap, p1, p2

115def _bandgap(e_skn, spin, direct):

116 """Helper function."""

117 ns, nk, nb = e_skn.shape

118 s1 = s2 = k1 = k2 = n1 = n2 = None

120 N_sk = (e_skn < 0.0).sum(2) # number of occupied bands

122 # Check for bands crossing the fermi-level

123 if ns == 1:

124 if N_sk[0].ptp() > 0:

125 return 0.0, (None, None, None), (None, None, None)

126 elif spin is None:

127 if (N_sk.ptp(axis=1) > 0).any():

128 return 0.0, (None, None, None), (None, None, None)

129 elif N_sk[spin].ptp() > 0:

130 return 0.0, (None, None, None), (None, None, None)

132 if (N_sk == 0).any() or (N_sk == nb).any():

133 raise ValueError('Too few bands!')

135 e_skn = np.array([[e_skn[s, k, N_sk[s, k] - 1:N_sk[s, k] + 1]

136 for k in range(nk)]

137 for s in range(ns)])

138 ev_sk = e_skn[:, :, 0] # valence band

139 ec_sk = e_skn[:, :, 1] # conduction band

141 if ns == 1:

142 s1 = 0

143 s2 = 0

144 gap, k1, k2 = find_gap(ev_sk[0], ec_sk[0], direct)

145 n1 = N_sk[0, 0] - 1

146 n2 = n1 + 1

147 return gap, (0, k1, n1), (0, k2, n2)

149 if spin is None:

150 gap, k1, k2 = find_gap(ev_sk.ravel(), ec_sk.ravel(), direct)

151 if direct:

152 # Check also spin flips:

153 for s in [0, 1]:

154 g, k, _ = find_gap(ev_sk[s], ec_sk[1 - s], direct)

155 if g < gap:

156 gap = g

157 k1 = k + nk * s

158 k2 = k + nk * (1 - s)

160 if gap > 0.0:

161 s1, k1 = divmod(k1, nk)

162 s2, k2 = divmod(k2, nk)

163 n1 = N_sk[s1, k1] - 1

164 n2 = N_sk[s2, k2]

165 return gap, (s1, k1, n1), (s2, k2, n2)

166 return 0.0, (None, None, None), (None, None, None)

168 gap, k1, k2 = find_gap(ev_sk[spin], ec_sk[spin], direct)

169 s1 = spin

170 s2 = spin

171 n1 = N_sk[s1, k1] - 1

172 n2 = n1 + 1

173 return gap, (s1, k1, n1), (s2, k2, n2)

176def find_gap(ev_k, ec_k, direct):

177 """Helper function."""

178 if direct:

179 gap_k = ec_k - ev_k

180 k = gap_k.argmin()

181 return gap_k[k], k, k

182 kv = ev_k.argmax()

183 kc = ec_k.argmin()

184 return ec_k[kc] - ev_k[kv], kv, kc