Coverage for /builds/ase/ase/ase/dft/dos.py: 79.41%

1# fmt: off 

2 

3from math import pi, sqrt 

4 

5import numpy as np 

6 

7from ase.dft.kpoints import get_monkhorst_pack_size_and_offset 

8from ase.parallel import world 

9from ase.utils.cext import cextension 

10 

11 

12class DOS: 

13 def __init__(self, calc, width=0.1, window=None, npts=401, comm=world): 

14 """Electronic Density Of States object. 

15 

16 calc: calculator object 

17 Any ASE compliant calculator object. 

18 width: float 

19 Width of guassian smearing. Use width=0.0 for linear tetrahedron 

20 interpolation. 

21 window: tuple of two float 

22 Use ``window=(emin, emax)``. If not specified, a window 

23 big enough to hold all the eigenvalues will be used. 

24 npts: int 

25 Number of points. 

26 comm: communicator object 

27 MPI communicator for lti_dos 

28 

29 """ 

30 

31 self.comm = comm 

32 self.npts = npts 

33 self.width = width 

34 self.w_k = calc.get_k_point_weights() 

35 self.nspins = calc.get_number_of_spins() 

36 self.e_skn = np.array([[calc.get_eigenvalues(kpt=k, spin=s) 

37 for k in range(len(self.w_k))] 

38 for s in range(self.nspins)]) 

39 try: # two Fermi levels 

40 for i, eF in enumerate(calc.get_fermi_level()): 

41 self.e_skn[i] -= eF 

42 except TypeError: # a single Fermi level 

43 self.e_skn -= calc.get_fermi_level() 

44 

45 if window is None: 

46 emin = None 

47 emax = None 

48 else: 

49 emin, emax = window 

50 

51 if emin is None: 

52 emin = self.e_skn.min() - 5 * self.width 

53 if emax is None: 

54 emax = self.e_skn.max() + 5 * self.width 

55 

56 self.energies = np.linspace(emin, emax, npts) 

57 

58 if width == 0.0: 

59 bzkpts = calc.get_bz_k_points() 

60 size, _offset = get_monkhorst_pack_size_and_offset(bzkpts) 

61 bz2ibz = calc.get_bz_to_ibz_map() 

62 shape = (self.nspins,) + tuple(size) + (-1,) 

63 self.e_skn = self.e_skn[:, bz2ibz].reshape(shape) 

64 self.cell = calc.atoms.cell 

65 

66 def get_energies(self): 

67 """Return the array of energies used to sample the DOS. 

68 

69 The energies are reported relative to the Fermi level. 

70 """ 

71 return self.energies 

72 

73 def delta(self, energy): 

74 """Return a delta-function centered at 'energy'.""" 

75 x = -((self.energies - energy) / self.width)**2 

76 return np.exp(x) / (sqrt(pi) * self.width) 

77 

78 def get_dos(self, spin=None): 

79 """Get array of DOS values. 

80 

81 The *spin* argument can be 0 or 1 (spin up or down) - if not 

82 specified, the total DOS is returned. 

83 """ 

84 

85 if spin is None: 

86 if self.nspins == 2: 

87 # Return the total DOS 

88 return self.get_dos(spin=0) + self.get_dos(spin=1) 

89 else: 

90 return 2 * self.get_dos(spin=0) 

91 elif spin == 1 and self.nspins == 1: 

92 # For an unpolarized calculation, spin up and down are equivalent 

93 spin = 0 

94 

95 if self.width == 0.0: 

96 dos = linear_tetrahedron_integration(self.cell, self.e_skn[spin], 

97 self.energies, comm=self.comm) 

98 return dos 

99 

100 dos = np.zeros(self.npts) 

101 for w, e_n in zip(self.w_k, self.e_skn[spin]): 

102 for e in e_n: 

103 dos += w * self.delta(e) 

104 return dos 

105 

106 

107def linear_tetrahedron_integration(cell, eigs, energies, 

108 weights=None, comm=world): 

109 """DOS from linear tetrahedron interpolation. 

110 

111 cell: 3x3 ndarray-like 

112 Unit cell. 

113 eigs: (n1, n2, n3, nbands)-shaped ndarray 

114 Eigenvalues on a Monkhorst-Pack grid (not reduced). 

115 energies: 1-d array-like 

116 Energies where the DOS is calculated (must be a uniform grid). 

117 weights: ndarray of shape (n1, n2, n3, nbands) or (n1, n2, n3, nbands, nw) 

118 Weights. Defaults to a (n1, n2, n3, nbands)-shaped ndarray 

119 filled with ones. Can also have an extra dimednsion if there are 

120 nw weights. 

121 comm: communicator object 

122 MPI communicator for lti_dos 

123 

124 Returns: 

125 

126 DOS as an ndarray of same length as energies or as an 

127 ndarray of shape (nw, len(energies)). 

128 

129 See: 

130 

131 Extensions of the tetrahedron method for evaluating 

132 spectral properties of solids, 

133 A. H. MacDonald, S. H. Vosko and P. T. Coleridge, 

134 1979 J. Phys. C: Solid State Phys. 12 2991, 

135 :doi:`10.1088/0022-3719/12/15/008` 

136 """ 

137 

138 from scipy.spatial import Delaunay 

139 

140 # Find the 6 tetrahedra: 

141 size = eigs.shape[:3] 

142 B = (np.linalg.inv(cell) / size).T 

143 indices = np.array([[i, j, k] 

144 for i in [0, 1] for j in [0, 1] for k in [0, 1]]) 

145 dt = Delaunay(np.dot(indices, B)) 

146 

147 if weights is None: 

148 weights = np.ones_like(eigs) 

149 

150 if weights.ndim == 4: 

151 extra_dimension_added = True 

152 weights = weights[:, :, :, :, np.newaxis] 

153 else: 

154 extra_dimension_added = False 

155 

156 nweights = weights.shape[4] 

157 dos = np.empty((nweights, len(energies))) 

158 

159 lti_dos(indices[dt.simplices], eigs, weights, energies, dos, comm) 

160 

161 dos /= np.prod(size) 

162 

163 if extra_dimension_added: 

164 return dos[0] 

165 return dos 

166 

167 

168@cextension 

169def lti_dos(simplices, eigs, weights, energies, dos, world): 

170 shape = eigs.shape[:3] 

171 nweights = weights.shape[-1] 

172 dos[:] = 0.0 

173 n = -1 

174 for index in np.indices(shape).reshape((3, -1)).T: 

175 n += 1 

176 if n % world.size != world.rank: 

177 continue 

178 i = ((index + simplices) % shape).T 

179 E = eigs[i[0], i[1], i[2]].reshape((4, -1)) 

180 W = weights[i[0], i[1], i[2]].reshape((4, -1, nweights)) 

181 for e, w in zip(E.T, W.transpose((1, 0, 2))): 

182 lti_dos1(e, w, energies, dos) 

183 

184 dos /= 6.0 

185 world.sum(dos) 

186 

187 

188def lti_dos1(e, w, energies, dos): 

189 i = e.argsort() 

190 e0, e1, e2, e3 = en = e[i] 

191 w = w[i] 

192 

193 zero = energies[0] 

194 if len(energies) > 1: 

195 de = energies[1] - zero 

196 nn = (np.floor((en - zero) / de).astype(int) + 1).clip(0, 

197 len(energies)) 

198 else: 

199 nn = (en > zero).astype(int) 

200 

201 n0, n1, n2, n3 = nn 

202 

203 if n1 > n0: 

204 s = slice(n0, n1) 

205 x = energies[s] - e0 

206 f10 = x / (e1 - e0) 

207 f20 = x / (e2 - e0) 

208 f30 = x / (e3 - e0) 

209 f01 = 1 - f10 

210 f02 = 1 - f20 

211 f03 = 1 - f30 

212 g = f20 * f30 / (e1 - e0) 

213 dos[:, s] += w.T.dot([f01 + f02 + f03, 

214 f10, 

215 f20, 

216 f30]) * g 

217 if n2 > n1: 

218 delta = e3 - e0 

219 s = slice(n1, n2) 

220 x = energies[s] 

221 f20 = (x - e0) / (e2 - e0) 

222 f30 = (x - e0) / (e3 - e0) 

223 f21 = (x - e1) / (e2 - e1) 

224 f31 = (x - e1) / (e3 - e1) 

225 f02 = 1 - f20 

226 f03 = 1 - f30 

227 f12 = 1 - f21 

228 f13 = 1 - f31 

229 g = 3 / delta * (f12 * f20 + f21 * f13) 

230 dos[:, s] += w.T.dot([g * f03 / 3 + f02 * f20 * f12 / delta, 

231 g * f12 / 3 + f13 * f13 * f21 / delta, 

232 g * f21 / 3 + f20 * f20 * f12 / delta, 

233 g * f30 / 3 + f31 * f13 * f21 / delta]) 

234 if n3 > n2: 

235 s = slice(n2, n3) 

236 x = energies[s] - e3 

237 f03 = x / (e0 - e3) 

238 f13 = x / (e1 - e3) 

239 f23 = x / (e2 - e3) 

240 f30 = 1 - f03 

241 f31 = 1 - f13 

242 f32 = 1 - f23 

243 g = f03 * f13 / (e3 - e2) 

244 dos[:, s] += w.T.dot([f03, 

245 f13, 

246 f23, 

247 f30 + f31 + f32]) * g