Coverage for /builds/hweiske/ase/ase/dft/bandgap.py: 85.71%

1import warnings 

2 

3import numpy as np 

4 

5 

6def get_band_gap(calc, direct=False, spin=None): 

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

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

9 ns = calc.get_number_of_spins() 

10 if ns == 2 and spin is None: 

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

12 return gap, k1, k2 

13 

14 

15def bandgap(calc=None, direct=False, spin=None, eigenvalues=None, efermi=None, 

16 output=None, kpts=None): 

17 """Calculates the band-gap. 

18 

19 Parameters: 

20 

21 calc: Calculator object 

22 Electronic structure calculator object. 

23 direct: bool 

24 Calculate direct band-gap. 

25 spin: int or None 

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

27 at a single spin-channel. 

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

29 Eigenvalues. 

30 efermi: float 

31 Fermi level (defaults to 0.0). 

32 

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

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

35 

36 Example: 

37 

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

39 Gap: 1.2 eV 

40 Transition (v -> c): 

41 [0.000, 0.000, 0.000] -> [0.500, 0.500, 0.000] 

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

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

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

45 Direct gap: 3.4 eV 

46 Transition at: [0.000, 0.000, 0.000] 

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

48 3.4 (0, 0, 3), (0, 0, 4) 

49 """ 

50 

51 if calc: 

52 kpts = calc.get_ibz_k_points() 

53 nk = len(kpts) 

54 ns = calc.get_number_of_spins() 

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

56 for k in range(nk)] 

57 for s in range(ns)]) 

58 if efermi is None: 

59 efermi = calc.get_fermi_level() 

60 

61 efermi = efermi or 0.0 

62 

63 e_skn = eigenvalues - efermi 

64 if eigenvalues.ndim == 2: 

65 e_skn = e_skn[np.newaxis] # spinors 

66 

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

68 raise ValueError('Bad eigenvalues!') 

69 

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

71 

72 if eigenvalues.ndim != 3: 

73 p1 = (k1, n1) 

74 p2 = (k2, n2) 

75 else: 

76 p1 = (s1, k1, n1) 

77 p2 = (s2, k2, n2) 

78 

79 return gap, p1, p2 

80 

81 

82def _bandgap(e_skn, spin, direct): 

83 """Helper function.""" 

84 ns, nk, nb = e_skn.shape 

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

86 

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

88 

89 # Check for bands crossing the fermi-level 

90 if ns == 1: 

91 if np.ptp(N_sk[0]) > 0: 

92 return 0.0, (None, None, None), (None, None, None) 

93 elif spin is None: 

94 if (np.ptp(N_sk, axis=1) > 0).any(): 

95 return 0.0, (None, None, None), (None, None, None) 

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

97 return 0.0, (None, None, None), (None, None, None) 

98 

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

100 raise ValueError('Too few bands!') 

101 

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

103 for k in range(nk)] 

104 for s in range(ns)]) 

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

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

107 

108 if ns == 1: 

109 s1 = 0 

110 s2 = 0 

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

112 n1 = N_sk[0, 0] - 1 

113 n2 = n1 + 1 

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

115 

116 if spin is None: 

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

118 if direct: 

119 # Check also spin flips: 

120 for s in [0, 1]: 

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

122 if g < gap: 

123 gap = g 

124 k1 = k + nk * s 

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

126 

127 if gap > 0.0: 

128 s1, k1 = divmod(k1, nk) 

129 s2, k2 = divmod(k2, nk) 

130 n1 = N_sk[s1, k1] - 1 

131 n2 = N_sk[s2, k2] 

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

133 return 0.0, (None, None, None), (None, None, None) 

134 

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

136 s1 = spin 

137 s2 = spin 

138 n1 = N_sk[s1, k1] - 1 

139 n2 = n1 + 1 

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

141 

142 

143def find_gap(ev_k, ec_k, direct): 

144 """Helper function.""" 

145 if direct: 

146 gap_k = ec_k - ev_k 

147 k = gap_k.argmin() 

148 return gap_k[k], k, k 

149 kv = ev_k.argmax() 

150 kc = ec_k.argmin() 

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