Coverage for /builds/kinetik161/ase/ase/dft/bandgap.py: 87.62%

1import functools 

2import warnings 

3 

4import numpy as np 

5 

6from ase.utils import IOContext 

7 

8 

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

10 warnings.warn('Please use ase.dft.bandgap.bandgap() instead!') 

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 

16 

17 

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

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

20 """Calculates the band-gap. 

21 

22 Parameters: 

23 

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. 

39 

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). 

42 

43 Example: 

44 

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 """ 

57 

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() 

67 

68 efermi = efermi or 0.0 

69 

70 e_skn = eigenvalues - efermi 

71 if eigenvalues.ndim == 2: 

72 e_skn = e_skn[np.newaxis] # spinors 

73 

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

75 raise ValueError('Bad eigenvalues!') 

76 

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

78 

79 with IOContext() as iocontext: 

80 fd = iocontext.openfile(output) 

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

82 

83 def skn(s, k, n): 

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

85 if kpts is None: 

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

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

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

89 

90 if spin is not None: 

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

92 if gap == 0.0: 

93 p('No gap') 

94 elif direct: 

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

96 if s1 == s2: 

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

98 else: 

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

100 else: 

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

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

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

104 

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) 

111 

112 return gap, p1, p2 

113 

114 

115def _bandgap(e_skn, spin, direct): 

116 """Helper function.""" 

117 ns, nk, nb = e_skn.shape 

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

119 

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

121 

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) 

131 

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

133 raise ValueError('Too few bands!') 

134 

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 

140 

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) 

148 

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) 

159 

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) 

167 

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) 

174 

175 

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