Coverage for /builds/kinetik161/ase/ase/dft/dos.py: 77.94%

1from math import pi, sqrt 

2 

3import numpy as np 

4 

5from ase.dft.kpoints import get_monkhorst_pack_size_and_offset 

6from ase.parallel import world 

7from ase.utils.cext import cextension 

8 

9 

10class DOS: 

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

12 """Electronic Density Of States object. 

13 

14 calc: calculator object 

15 Any ASE compliant calculator object. 

16 width: float 

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

18 interpolation. 

19 window: tuple of two float 

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

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

22 npts: int 

23 Number of points. 

24 comm: communicator object 

25 MPI communicator for lti_dos 

26 

27 """ 

28 

29 self.comm = comm 

30 self.npts = npts 

31 self.width = width 

32 self.w_k = calc.get_k_point_weights() 

33 self.nspins = calc.get_number_of_spins() 

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

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

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

37 try: # two Fermi levels 

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

39 self.e_skn[i] -= eF 

40 except TypeError: # a single Fermi level 

41 self.e_skn -= calc.get_fermi_level() 

42 

43 if window is None: 

44 emin = None 

45 emax = None 

46 else: 

47 emin, emax = window 

48 

49 if emin is None: 

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

51 if emax is None: 

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

53 

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

55 

56 if width == 0.0: 

57 bzkpts = calc.get_bz_k_points() 

58 size, offset = get_monkhorst_pack_size_and_offset(bzkpts) 

59 bz2ibz = calc.get_bz_to_ibz_map() 

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

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

62 self.cell = calc.atoms.cell 

63 

64 def get_energies(self): 

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

66 

67 The energies are reported relative to the Fermi level. 

68 """ 

69 return self.energies 

70 

71 def delta(self, energy): 

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

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

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

75 

76 def get_dos(self, spin=None): 

77 """Get array of DOS values. 

78 

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

80 specified, the total DOS is returned. 

81 """ 

82 

83 if spin is None: 

84 if self.nspins == 2: 

85 # Return the total DOS 

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

87 else: 

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

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

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

91 spin = 0 

92 

93 if self.width == 0.0: 

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

95 self.energies, comm=self.comm) 

96 return dos 

97 

98 dos = np.zeros(self.npts) 

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

100 for e in e_n: 

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

102 return dos 

103 

104 

105def linear_tetrahedron_integration(cell, eigs, energies, 

106 weights=None, comm=world): 

107 """DOS from linear tetrahedron interpolation. 

108 

109 cell: 3x3 ndarray-like 

110 Unit cell. 

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

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

113 energies: 1-d array-like 

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

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

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

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

118 nw weights. 

119 comm: communicator object 

120 MPI communicator for lti_dos 

121 

122 Returns: 

123 

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

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

126 

127 See: 

128 

129 Extensions of the tetrahedron method for evaluating 

130 spectral properties of solids, 

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

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

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

134 """ 

135 

136 from scipy.spatial import Delaunay 

137 

138 # Find the 6 tetrahedra: 

139 size = eigs.shape[:3] 

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

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

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

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

144 

145 if weights is None: 

146 weights = np.ones_like(eigs) 

147 

148 if weights.ndim == 4: 

149 extra_dimension_added = True 

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

151 else: 

152 extra_dimension_added = False 

153 

154 nweights = weights.shape[4] 

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

156 

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

158 

159 dos /= np.prod(size) 

160 

161 if extra_dimension_added: 

162 return dos[0] 

163 return dos 

164 

165 

166@cextension 

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

168 shape = eigs.shape[:3] 

169 nweights = weights.shape[-1] 

170 dos[:] = 0.0 

171 n = -1 

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

173 n += 1 

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

175 continue 

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

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

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

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

180 lti_dos1(e, w, energies, dos) 

181 

182 dos /= 6.0 

183 world.sum(dos) 

184 

185 

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

187 i = e.argsort() 

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

189 w = w[i] 

190 

191 zero = energies[0] 

192 if len(energies) > 1: 

193 de = energies[1] - zero 

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

195 len(energies)) 

196 else: 

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

198 

199 n0, n1, n2, n3 = nn 

200 

201 if n1 > n0: 

202 s = slice(n0, n1) 

203 x = energies[s] - e0 

204 f10 = x / (e1 - e0) 

205 f20 = x / (e2 - e0) 

206 f30 = x / (e3 - e0) 

207 f01 = 1 - f10 

208 f02 = 1 - f20 

209 f03 = 1 - f30 

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

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

212 f10, 

213 f20, 

214 f30]) * g 

215 if n2 > n1: 

216 delta = e3 - e0 

217 s = slice(n1, n2) 

218 x = energies[s] 

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

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

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

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

223 f02 = 1 - f20 

224 f03 = 1 - f30 

225 f12 = 1 - f21 

226 f13 = 1 - f31 

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

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

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

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

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

232 if n3 > n2: 

233 s = slice(n2, n3) 

234 x = energies[s] - e3 

235 f03 = x / (e0 - e3) 

236 f13 = x / (e1 - e3) 

237 f23 = x / (e2 - e3) 

238 f30 = 1 - f03 

239 f31 = 1 - f13 

240 f32 = 1 - f23 

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

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

243 f13, 

244 f23, 

245 f30 + f31 + f32]) * g