Coverage for eminus/occupations.py: 98.28%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Determine occupations for atomic systems from simple inputs.""" 

4 

5import dataclasses 

6import numbers 

7 

8import numpy as np 

9 

10from .logger import log 

11from .tools import fermi_distribution, get_Efermi 

12 

13 

14@dataclasses.dataclass 

15class Occupations: 

16 """Occupations class to save electronic state information in one place. 

17 

18 The attribute Nelec has to be given first after instantiation. 

19 """ 

20 

21 # Set the private variable for the attributes that are properties. 

22 _Nelec: int = 0 #: Number of electrons. 

23 _Nspin: int = 0 #: Number of spin states. 

24 _spin: float = 0 #: Number of unpaired electrons. 

25 _charge: int = 0 #: System charge. 

26 _Nstate: int = 0 #: Number of states. 

27 _Nempty: int = 0 #: Number of empty states. 

28 _Nk: int = 1 #: Number of k-points. 

29 _bands: int = 0 #: Number of bands. 

30 _smearing: float = 0 #: Smearing width in Hartree. 

31 is_filled: bool = False #: Determines the Occupations object fill status. 

32 

33 # ### Class properties ### 

34 

35 @property 

36 def Nelec(self): 

37 """Number of electrons.""" 

38 return self._Nelec 

39 

40 @Nelec.setter 

41 def Nelec(self, value): 

42 # Only update Nelec if it actually gets updated 

43 if self._Nelec != int(value): 

44 self._Nelec = int(value) 

45 self.is_filled = False 

46 

47 @property 

48 def Nspin(self): 

49 """Number of spin states.""" 

50 return self._Nspin 

51 

52 @Nspin.setter 

53 def Nspin(self, value): 

54 # Use a spin-paired calculation for an even number of electrons 

55 if value is None: 

56 if self.Nelec % 2 == 0 and self.spin == 0: 

57 value = 1 

58 else: 

59 value = 2 

60 # Only update if needed 

61 if self._Nspin != int(value): 

62 self._Nspin = int(value) 

63 self.spin = self.spin 

64 self.is_filled = False 

65 if hasattr(self, "f"): 

66 log.warning("Reset previously set fillings.") 

67 

68 @property 

69 def spin(self): 

70 """Number of unpaired electrons.""" 

71 return self._spin 

72 

73 @spin.setter 

74 def spin(self, value): 

75 # If no spin is given try to determine it from the number of electrons 

76 if value is None: 

77 if self.Nelec % 2 == 0: 

78 value = 0 

79 else: 

80 value = 1 

81 if self._spin != value: 

82 self.is_filled = False 

83 # We have no spin in the spin-paired case 

84 if self.Nspin == 1: 

85 self._spin = 0 

86 else: 

87 self._spin = value 

88 

89 @property 

90 def charge(self): 

91 """System charge.""" 

92 return self._charge 

93 

94 @charge.setter 

95 def charge(self, value): 

96 # If we set a charge via this setter update the number of electrons 

97 if hasattr(self, "Nelec"): 

98 self.Nelec += self._charge - int(value) 

99 if self._charge != int(value): 

100 self._charge = int(value) 

101 self.is_filled = False 

102 

103 @property 

104 def f(self): 

105 """Occupation numbers per state.""" 

106 return self._f 

107 

108 @f.setter 

109 def f(self, value): 

110 # Make sure the occupations are in a two-dimensional array 

111 if isinstance(value, (list, tuple, np.ndarray)): 

112 value = np.atleast_2d(value) 

113 self.is_filled = False 

114 # This setter will only be called when explicitly setting f 

115 # Call the fill function in that case 

116 self.fill(value) 

117 

118 @property 

119 def Nk(self): 

120 """Number of k-points.""" 

121 return self._Nk 

122 

123 @Nk.setter 

124 def Nk(self, value): 

125 self._Nk = int(value) 

126 self.is_filled = False 

127 

128 @property 

129 def wk(self): 

130 """k-point weights.""" 

131 return self._wk 

132 

133 @wk.setter 

134 def wk(self, value): 

135 self._wk = np.asarray(value) 

136 self._Nk = len(self._wk) 

137 self.is_filled = False 

138 

139 @property 

140 def bands(self): 

141 """Total number of bands.""" 

142 return self._bands 

143 

144 @bands.setter 

145 def bands(self, value): 

146 self._bands = int(value) 

147 self.is_filled = False 

148 

149 @property 

150 def smearing(self): 

151 """Smearing width in Hartree.""" 

152 return self._smearing 

153 

154 @smearing.setter 

155 def smearing(self, value): 

156 if value < 0: 

157 log.error("The smearing width can not be negative.") 

158 self._smearing = value 

159 self.is_filled = False 

160 if self.Nempty > 0: 

161 log.warning("Empty states with smearing enabled found.") 

162 

163 @property 

164 def magnetization(self): 

165 """Magnetization from occupation numbers.""" 

166 # There is no magnetization in the spin-paired case 

167 if self.Nspin == 1: 

168 return 0 

169 if not self.is_filled: 

170 log.warning("Can not calculate magnetization for unfilled occupations.") 

171 return 0 

172 return np.sum(self.wk * (self.f[:, 0] - self.f[:, 1])) / np.sum(self.wk * self.f) 

173 

174 @magnetization.setter 

175 def magnetization(self, value): 

176 if value is not None: 

177 if self.is_filled: 

178 log.warning("Reset previously set fillings.") 

179 self.is_filled = False 

180 self.fill(None, value) 

181 

182 # ### Read-only properties ### 

183 

184 @property 

185 def multiplicity(self): 

186 """Multiplicity, i.e., 2S+1.""" 

187 return self.spin + 1 

188 

189 @property 

190 def Nstate(self): 

191 """Number of states.""" 

192 return self._Nstate 

193 

194 @property 

195 def Nempty(self): 

196 """Number of empty states.""" 

197 return self._Nempty 

198 

199 @property 

200 def F(self): 

201 """Diagonal matrices of f per k-point and spin.""" 

202 return [[np.diag(f) for f in f_spin] for f_spin in self.f] 

203 

204 # ### Class methods ### 

205 

206 def fill(self, f=None, magnetization=None): 

207 """Fill the states of the object. 

208 

209 Keyword Args: 

210 f: Fillings. 

211 magnetization: Magnetization. 

212 """ 

213 # Do nothing if the object is already filled 

214 if self.is_filled: 

215 return self 

216 # If no f is given just use the standard fillings: 2 for restricted and 1 for unrestricted 

217 if f is None: 

218 f = 2 / self.Nspin 

219 self._update_from_fillings(f, magnetization) 

220 self.is_filled = True 

221 return self 

222 

223 kernel = fill 

224 

225 def _update_from_fillings(self, value, magnetization): 

226 """Update fillings. 

227 

228 Args: 

229 value: Fillings. 

230 magnetization: Magnetization. 

231 """ 

232 # Do not use the setter methods in this place to not trigger the setter effects 

233 # If f is a number use this occupation for all states 

234 if isinstance(value, numbers.Real) or isinstance(magnetization, numbers.Real): 

235 # Do not leave the states array empty when no electrons are present 

236 if self.Nelec <= 0: 

237 self._Nstate = 1 

238 self._f = np.zeros((self.Nk, self.Nspin, 1)) 

239 # Always use the fractional fillings method if a magnetization is given 

240 elif isinstance(magnetization, numbers.Real): 

241 self._fractional_fillings(value, magnetization) 

242 elif self.Nspin == 1 or self.Nelec % 2 == self.spin % 2: 

243 self._integer_fillings(value) 

244 else: 

245 log.warning("Using fractional occupations.") 

246 self._fractional_fillings(value) 

247 # If f is an array redetermine all attributes from it 

248 else: 

249 self._f = value 

250 self._Nspin = len(value) 

251 self._Nstate = len(value[0]) 

252 self._charge += self.Nelec - int(np.sum(value)) 

253 self.Nelec = int(np.sum(value)) 

254 if self.Nspin == 1: 

255 self.spin = 0 

256 else: 

257 self.spin = abs(np.sum(self.f[0] - self.f[1])) 

258 

259 def _integer_fillings(self, f): 

260 """Update fillings while maintaining integer occupations numbers. 

261 

262 Args: 

263 f: Fillings. 

264 """ 

265 # Determine the electrons per spin channel 

266 Nup = self.Nelec // self.Nspin + self.spin // self.Nspin + self.Nelec % self.Nspin 

267 Ndw = self.Nelec // self.Nspin - self.spin // self.Nspin 

268 elecs = np.array([Nup, Ndw]) 

269 

270 # Get the number of states 

271 Nstate = int(np.ceil(max(elecs / f))) 

272 # If no bands are set, set them now 

273 if self.bands == 0: 

274 self.bands = Nstate 

275 self._Nstate = Nstate 

276 if self.bands < Nstate: 

277 log.error("Number of bands is smaller than the number of valence electrons.") 

278 if self.bands == Nstate and self.smearing > 0: 

279 log.warning("Smearing has been enabled but no extra bands have been set.") 

280 # Set the number of empty states if no smearing is used 

281 if self.smearing == 0: 

282 self._Nstate = Nstate 

283 self._Nempty = self.bands - Nstate 

284 

285 # Simply build the occupations array 

286 self._f = f * np.ones((self.Nspin, Nstate), dtype=int) 

287 # If we have filled too much correct it in the second spin channel 

288 # The following procedure ensures that we have no negative fillings 

289 rest = np.sum(self._f) - self.Nelec 

290 i = 1 

291 while rest > 0: 

292 # If the occupation is greater than the rest we are done with removing electrons 

293 if self._f[-1, -i] >= rest: 

294 self._f[-1, -i] -= rest 

295 break 

296 # Otherwise zero out the state, update the rest, and move to the next state 

297 rest -= self._f[-1, -i] 

298 self._f[-1, -i] = 0 

299 i += 1 

300 

301 if self.smearing > 0: 

302 # Append extra states 

303 self._f = np.hstack((self._f, np.zeros((self.Nspin, self.bands - Nstate)))) 

304 self._Nstate = self.bands 

305 self._f = np.vstack([[self._f]] * self.Nk) 

306 

307 def _fractional_fillings(self, f, magnetization=None): 

308 """Update fillings while allowing fractional occupation numbers. 

309 

310 Args: 

311 f: Fillings. 

312 

313 Keyword Args: 

314 magnetization: Magnetization. 

315 """ 

316 # Determine the electrons per spin channel 

317 if magnetization is None: 

318 Nup = self.Nelec / self.Nspin + self.spin / self.Nspin 

319 Ndw = self.Nelec / self.Nspin - self.spin / self.Nspin 

320 else: 

321 # M = (Nup - Ndw) / N = (Nup - N + Nup) / N 

322 Nup = (magnetization * self.Nelec + self.Nelec) / 2 

323 Ndw = self.Nelec - Nup 

324 self.spin = abs(Nup - Ndw) 

325 elecs = np.array([Nup, Ndw]) 

326 

327 # Get the number of states 

328 Nstate = int(np.ceil(max(elecs / f))) 

329 # If no bands are set, set them now 

330 if self.bands == 0: 

331 self.bands = Nstate 

332 self._Nstate = Nstate 

333 if self.bands < Nstate: 

334 log.error("Number of bands is smaller than the number of valence electrons.") 

335 if self.bands == Nstate and self.smearing > 0: 

336 log.warning("Smearing has been enabled but no extra bands have been set.") 

337 # Set the number of empty states if no smearing is used 

338 if self.smearing == 0: 

339 self._Nstate = Nstate 

340 self._Nempty = self.bands - Nstate 

341 

342 # Simply build the occupations array 

343 self._f = f * np.ones((self.Nspin, Nstate)) 

344 # If we have filled too much correct it in both spin channels 

345 # The following procedure ensures that we have no negative fillings 

346 rest = np.sum(self._f, axis=1) - elecs 

347 for spin in range(self.Nspin): 

348 i = 1 

349 while rest[spin] > 0: 

350 # If the occupation is greater than the rest we are done with removing electrons 

351 if self._f[spin, -i] >= rest[spin]: 

352 self._f[spin, -i] -= rest[spin] 

353 break 

354 # Otherwise zero out the state, update the rest, and move to the next state 

355 rest[spin] -= self._f[spin, -i] 

356 self._f[spin, -i] = 0 

357 i += 1 

358 

359 if self.smearing > 0: 

360 # Append extra states 

361 self._f = np.hstack((self._f, np.zeros((self.Nspin, self.bands - Nstate)))) 

362 self._Nstate = self.bands 

363 self._f = np.vstack([[self._f]] * self.Nk) 

364 

365 def __repr__(self): 

366 """Print the parameters stored in the Occupations object.""" 

367 return ( 

368 f"Spin handling: {'un' if self.Nspin == 2 else ''}restricted\n" 

369 f"Number of electrons: {self.Nelec}\n" 

370 f"Spin: {self.spin}\n" 

371 f"Charge: {self.charge}\n" 

372 f"Number of bands: {self.bands}\n" 

373 f"Number of states: {self.Nstate}\n" 

374 f"Number of empty states: {self.Nempty}\n" 

375 f"Number of k-points: {self.Nk}\n" 

376 f"Smearing width: {self.smearing} Eh\n" 

377 f"Fillings: \n{self.f if self.is_filled else 'Not filled'}" 

378 ) 

379 

380 def smear(self, epsilon): 

381 """Update fillings according to a Fermi distribution. 

382 

383 Args: 

384 epsilon: Eigenenergies. 

385 

386 Returns: 

387 Efermi: Fermi energy. 

388 """ 

389 if self.smearing == 0: 

390 log.info("Smearing is set to zero, nothing to do.") 

391 return 0 

392 

393 Efermi = get_Efermi(self, epsilon) 

394 self._f = fermi_distribution(epsilon, Efermi, self.smearing) * 2 / self.Nspin 

395 return Efermi