1"""Smearing functions and occupation number calculators.""" 2 3import warnings 4from math import pi, nan, inf 5from typing import List, Tuple, NamedTuple, Any, Callable, Dict, cast 6import numpy as np 7from scipy.special import erf 8from ase.units import Ha 9 10from gpaw.band_descriptor import BandDescriptor 11from gpaw.mpi import serial_comm, broadcast_float 12from gpaw.typing import Array1D, Array2D 13 14# typehints: 15MPICommunicator = Any 16 17 18class ParallelLayout(NamedTuple): 19 """Collection of parallel stuff.""" 20 bd: BandDescriptor 21 kpt_comm: MPICommunicator 22 domain_comm: MPICommunicator 23 24 25def fermi_dirac(eig: np.ndarray, 26 fermi_level: float, 27 width: float) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: 28 """Fermi-Dirac distribution function. 29 30 >>> f, _, _ = fermi_dirac(0.0, 0.0, 0.1) 31 >>> f 32 0.5 33 """ 34 x = (eig - fermi_level) / width 35 x = np.clip(x, -100, 100) 36 y = np.exp(x) 37 z = y + 1.0 38 f = 1.0 / z 39 dfde = (f - f**2) / width 40 y *= x 41 y /= z 42 y -= np.log(z) 43 e_entropy = y * width 44 return f, dfde, e_entropy 45 46 47def marzari_vanderbilt(eig: np.ndarray, 48 fermi_level: float, 49 width: float) -> Tuple[np.ndarray, 50 np.ndarray, 51 np.ndarray]: 52 """Marzari-Vanderbilt distribution (cold smearing). 53 54 See: https://doi.org/10.1103/PhysRevLett.82.3296 55 """ 56 x = (eig - fermi_level) / width 57 expterm = np.exp(-(x + (1 / np.sqrt(2)))**2) 58 f = expterm / np.sqrt(2 * np.pi) + 0.5 * (1 - erf(1. / np.sqrt(2) + x)) 59 dfde = expterm * (2 + np.sqrt(2) * x) / np.sqrt(np.pi) / width 60 s = expterm * (1 + np.sqrt(2) * x) / (2 * np.sqrt(np.pi)) 61 e_entropy = -s * width 62 return f, dfde, e_entropy 63 64 65def methfessel_paxton(eig: np.ndarray, 66 fermi_level: float, 67 width: float, 68 order: int = 0) -> Tuple[np.ndarray, 69 np.ndarray, 70 np.ndarray]: 71 """Methfessel-Paxton distribution.""" 72 x = (eig - fermi_level) / width 73 f = 0.5 * (1 - erf(x)) 74 for i in range(order): 75 f += (coff_function(i + 1) * 76 hermite_poly(2 * i + 1, x) * np.exp(-x**2)) 77 dfde = 1 / np.sqrt(pi) * np.exp(-x**2) 78 for i in range(order): 79 dfde += (coff_function(i + 1) * 80 hermite_poly(2 * i + 2, x) * np.exp(-x**2)) 81 dfde *= 1.0 / width 82 e_entropy = (0.5 * coff_function(order) * 83 hermite_poly(2 * order, x) * np.exp(-x**2)) 84 e_entropy = -e_entropy * width 85 return f, dfde, e_entropy 86 87 88def coff_function(n): 89 return (-1)**n / (np.product(np.arange(1, n + 1)) * 90 4**n * np.sqrt(np.pi)) 91 92 93def hermite_poly(n, x): 94 if n == 0: 95 return 1 96 elif n == 1: 97 return 2 * x 98 else: 99 return (2 * x * hermite_poly(n - 1, x) - 100 2 * (n - 1) * hermite_poly(n - 2, x)) 101 102 103class OccupationNumberCalculator: 104 """Base class for all occupation number calculators.""" 105 name = 'unknown' 106 extrapolate_factor: float 107 108 def __init__(self, 109 parallel_layout: ParallelLayout = None): 110 """Object for calculating fermi level(s) and occupation numbers. 111 112 If fixmagmom=True then the fixed_magmom_value attribute must be set 113 and two fermi levels will be calculated. 114 """ 115 if parallel_layout is None: 116 parallel_layout = ParallelLayout(BandDescriptor(1), 117 serial_comm, 118 serial_comm) 119 self.bd = parallel_layout.bd 120 self.kpt_comm = parallel_layout.kpt_comm 121 self.domain_comm = parallel_layout.domain_comm 122 123 @property 124 def parallel_layout(self) -> ParallelLayout: 125 return ParallelLayout(self.bd, self.kpt_comm, self.domain_comm) 126 127 def todict(self): 128 return {'name': self.name} 129 130 def copy(self, 131 parallel_layout: ParallelLayout = None, 132 bz2ibzmap: List[int] = None 133 ) -> 'OccupationNumberCalculator': 134 return create_occ_calc( 135 self.todict(), 136 parallel_layout=parallel_layout or self.parallel_layout) 137 138 def calculate(self, 139 nelectrons: float, 140 eigenvalues: List[List[float]], 141 weights: List[float], 142 fermi_levels_guess: List[float] = None 143 ) -> Tuple[Array2D, 144 List[float], 145 float]: 146 """Calculate occupation numbers and fermi level(s) from eigenvalues. 147 148 nelectrons: 149 Number of electrons. 150 eigenvalues: ndarray, shape=(nibzkpts, nbands) 151 Eigenvalues in Hartree. 152 weights: ndarray, shape=(nibzkpts,) 153 Weights of k-points in IBZ (must sum to 1). 154 parallel: 155 Parallel distribution of eigenvalues. 156 fermi_level_guesses: 157 Optional guess(es) at fermi level(s). 158 159 Returns a tuple containing: 160 161 * occupation numbers (in the range 0 to 1) 162 * fermi-level in Hartree 163 * entropy as -S*T in Hartree 164 165 >>> occ = ZeroWidth() 166 >>> occ.calculate(1, [[0, 1]], [1]) 167 (array([[1., 0.]]), [0.5], 0.0) 168 """ 169 170 eig_qn = [np.asarray(eig_n) for eig_n in eigenvalues] 171 weight_q = np.asarray(weights) 172 173 if fermi_levels_guess is None: 174 fermi_levels_guess = [nan] 175 176 f_qn = np.empty((len(weight_q), len(eig_qn[0]))) 177 178 result = np.empty(2) 179 180 if self.domain_comm.rank == 0: 181 # Let the master domain do the work and broadcast results: 182 result[:] = self._calculate( 183 nelectrons, eig_qn, weight_q, f_qn, fermi_levels_guess[0]) 184 185 self.domain_comm.broadcast(result, 0) 186 self.domain_comm.broadcast(f_qn, 0) 187 188 fermi_level, e_entropy = result 189 return f_qn, [fermi_level], e_entropy 190 191 def _calculate(self, 192 nelectrons: float, 193 eig_qn: List[Array1D], 194 weight_q: Array1D, 195 f_qn: Array2D, 196 fermi_level_guess: float) -> Tuple[float, float]: 197 raise NotImplementedError 198 199 200class FixMagneticMomentOccupationNumberCalculator(OccupationNumberCalculator): 201 """Base class for all occupation number objects.""" 202 name = 'fixmagmom' 203 204 def __init__(self, 205 occ: OccupationNumberCalculator, 206 magmom: float): 207 """Object for calculating fermi level(s) and occupation numbers. 208 209 If fixmagmom=True then the fixed_magmom_value attribute must be set 210 and two fermi levels will be calculated. 211 """ 212 self.occ = occ 213 self.fixed_magmom_value = magmom 214 self.extrapolate_factor = occ.extrapolate_factor 215 216 def todict(self): 217 dct = self.occ.todict() 218 dct['fixmagmom'] = True 219 return dct 220 221 def calculate(self, 222 nelectrons: float, 223 eigenvalues: List[List[float]], 224 weights: List[float], 225 fermi_levels_guess: List[float] = None 226 ) -> Tuple[Array2D, 227 List[float], 228 float]: 229 230 magmom = self.fixed_magmom_value 231 232 if fermi_levels_guess is None: 233 fermi_levels_guess = [nan, nan] 234 235 f1_qn, fermi_levels1, e_entropy1 = self.occ.calculate( 236 (nelectrons + magmom) / 2, 237 eigenvalues[::2], 238 weights[::2], 239 fermi_levels_guess[:1]) 240 241 f2_qn, fermi_levels2, e_entropy2 = self.occ.calculate( 242 (nelectrons - magmom) / 2, 243 eigenvalues[1::2], 244 weights[1::2], 245 fermi_levels_guess[1:]) 246 247 f_qn = [] 248 for f1_n, f2_n in zip(f1_qn, f2_qn): 249 f_qn += [f1_n, f2_n] 250 251 return (np.array(f_qn), 252 fermi_levels1 + fermi_levels2, 253 e_entropy1 + e_entropy2) 254 255 256class SmoothDistribution(OccupationNumberCalculator): 257 """Base class for Fermi-Dirac and other smooth distributions.""" 258 def __init__(self, width: float, parallel_layout: ParallelLayout = None): 259 """Smooth distribution. 260 261 width: float 262 Width of distribution in eV. 263 fixmagmom: bool 264 Fix spin moment calculations. A separate Fermi level for 265 spin up and down electrons is found. 266 """ 267 268 self._width = width 269 OccupationNumberCalculator.__init__(self, parallel_layout) 270 271 def todict(self): 272 return {'name': self.name, 'width': self._width} 273 274 def _calculate(self, 275 nelectrons, 276 eig_qn, 277 weight_q, 278 f_qn, 279 fermi_level_guess): 280 281 if np.isnan(fermi_level_guess) or self._width == 0.0: 282 zero = ZeroWidth(self.parallel_layout) 283 fermi_level_guess, _ = zero._calculate( 284 nelectrons, eig_qn, weight_q, f_qn) 285 if self._width == 0.0 or np.isinf(fermi_level_guess): 286 return fermi_level_guess, 0.0 287 288 x = fermi_level_guess 289 290 data = np.empty(3) 291 292 def func(x, data=data): 293 data[:] = 0.0 294 for eig_n, weight, f_n in zip(eig_qn, weight_q, f_qn): 295 f_n[:], dfde_n, e_entropy_n = self.distribution(eig_n, x) 296 data += [weight * x_n.sum() 297 for x_n in [f_n, dfde_n, e_entropy_n]] 298 self.bd.comm.sum(data) 299 self.kpt_comm.sum(data) 300 f, dfde = data[:2] 301 df = f - nelectrons 302 return df, dfde 303 304 fermi_level, niter = findroot(func, x) 305 306 e_entropy = data[2] 307 308 return fermi_level, e_entropy 309 310 311class FermiDiracCalculator(SmoothDistribution): 312 name = 'fermi-dirac' 313 extrapolate_factor = -0.5 314 315 def distribution(self, 316 eig_n: np.ndarray, 317 fermi_level: float) -> Tuple[np.ndarray, 318 np.ndarray, 319 np.ndarray]: 320 return fermi_dirac(eig_n, fermi_level, self._width) 321 322 def __str__(self): 323 return f'Fermi-Dirac: width={self._width:.4f} eV\n' 324 325 326class MarzariVanderbiltCalculator(SmoothDistribution): 327 name = 'marzari-vanderbilt' 328 # According to Nicola Marzari, one should not extrapolate M-V energies 329 # https://lists.quantum-espresso.org/pipermail/users/2005-October/003170.html 330 extrapolate_factor = 0.0 331 332 def distribution(self, eig_n, fermi_level): 333 return marzari_vanderbilt(eig_n, fermi_level, self._width) 334 335 def __str__(self): 336 return f'Marzari-Vanderbilt: width={self._width:.4f} eV\n' 337 338 339class MethfesselPaxtonCalculator(SmoothDistribution): 340 name = 'methfessel_paxton' 341 342 def __init__(self, width, order=0, parallel_layout: ParallelLayout = None): 343 SmoothDistribution.__init__(self, width, parallel_layout) 344 self.order = order 345 self.extrapolate_factor = -1.0 / (self.order + 2) 346 347 def todict(self): 348 dct = SmoothDistribution.todict(self) 349 dct['order'] = self.order 350 return dct 351 352 def __str__(self): 353 return (f'Methfessel-Paxton: width={self._width:.4f} eV, ' + 354 f'order={self.order}\n') 355 356 def distribution(self, eig_n, fermi_level): 357 return methfessel_paxton(eig_n, fermi_level, self._width, self.order) 358 359 360def findroot(func: Callable[[float], Tuple[float, float]], 361 x: float, 362 tol: float = 1e-10) -> Tuple[float, int]: 363 """Function used for locating Fermi level. 364 365 The function should return a (value, derivative) tuple: 366 367 >>> x, _ = findroot(lambda x: (x, 1.0), 1.0) 368 >>> assert abs(x) < 1e-10 369 """ 370 xmin = -np.inf 371 xmax = np.inf 372 373 # Try 10 step using the gradient: 374 niter = 0 375 while True: 376 f, dfdx = func(x) 377 if abs(f) < tol: 378 return x, niter 379 if f < 0.0 and x > xmin: 380 xmin = x 381 elif f > 0.0 and x < xmax: 382 xmax = x 383 dx = -f / max(dfdx, 1e-18) 384 if niter == 10 or abs(dx) > 0.3 or not (xmin < x + dx < xmax): 385 break # try bisection 386 x += dx 387 niter += 1 388 389 # Bracket the solution: 390 if not np.isfinite(xmin): 391 xmin = x 392 fmin = f 393 step = 0.01 394 while fmin > tol: 395 xmin -= step 396 fmin = func(xmin)[0] 397 step *= 2 398 399 if not np.isfinite(xmax): 400 xmax = x 401 fmax = f 402 step = 0.01 403 while fmax < 0: 404 xmax += step 405 fmax = func(xmax)[0] 406 step *= 2 407 408 # Bisect: 409 while True: 410 x = (xmin + xmax) / 2 411 f = func(x)[0] 412 if abs(f) < tol: 413 return x, niter 414 if f > 0: 415 xmax = x 416 else: 417 xmin = x 418 niter += 1 419 assert niter < 1000 420 421 422def collect_eigelvalues(eig_qn: np.ndarray, 423 weight_q: np.ndarray, 424 bd: BandDescriptor, 425 kpt_comm: MPICommunicator) -> Tuple[np.ndarray, 426 np.ndarray, 427 np.ndarray]: 428 """Collect eigenvalues from bd.comm and kpt_comm.""" 429 nkpts_r = np.zeros(kpt_comm.size, int) 430 nkpts_r[kpt_comm.rank] = len(weight_q) 431 kpt_comm.sum(nkpts_r) 432 nk = cast(int, nkpts_r.sum()) 433 weight_k = np.zeros(nk) 434 k1 = nkpts_r[:kpt_comm.rank].sum() 435 k2 = k1 + len(weight_q) 436 weight_k[k1:k2] = weight_q 437 kpt_comm.sum(weight_k, 0) 438 439 eig_kn: Array2D = np.zeros((0, 0)) 440 k = 0 441 for rank, nkpts in enumerate(nkpts_r): 442 for q in range(nkpts): 443 if rank == kpt_comm.rank: 444 eig_n = eig_qn[q] 445 eig_n = bd.collect(eig_n) 446 if bd.comm.rank == 0: 447 if kpt_comm.rank == 0: 448 if k == 0: 449 eig_kn = np.empty((nk, len(eig_n))) 450 if rank == 0: 451 eig_kn[k] = eig_n 452 else: 453 kpt_comm.receive(eig_kn[k], rank) 454 elif rank == kpt_comm.rank: 455 kpt_comm.send(eig_n, 0) 456 k += 1 457 return eig_kn, weight_k, nkpts_r 458 459 460def distribute_occupation_numbers(f_kn: np.ndarray, # input 461 f_qn: np.ndarray, # output 462 nkpts_r: np.ndarray, 463 bd: BandDescriptor, 464 kpt_comm: MPICommunicator) -> None: 465 """Distribute occupation numbers over bd.comm and kpt_comm.""" 466 k = 0 467 for rank, nkpts in enumerate(nkpts_r): 468 for q in range(nkpts): 469 if kpt_comm.rank == 0: 470 if rank == 0: 471 if bd.comm.size == 1: 472 f_qn[q] = f_kn[k] 473 else: 474 bd.distribute(None if f_kn is None else f_kn[k], 475 f_qn[q]) 476 elif f_kn is not None: 477 kpt_comm.send(f_kn[k], rank) 478 elif rank == kpt_comm.rank: 479 if bd.comm.size == 1: 480 kpt_comm.receive(f_qn[q], 0) 481 else: 482 if bd.comm.rank == 0: 483 f_n = bd.empty(global_array=True) 484 kpt_comm.receive(f_n, 0) 485 else: 486 f_n = None 487 bd.distribute(f_n, f_qn[q]) 488 k += 1 489 490 491class ZeroWidth(OccupationNumberCalculator): 492 name = 'zero-width' 493 extrapolate_factor = 0.0 494 495 def todict(self): 496 return {'width': 0.0} 497 498 def __str__(self): 499 return 'width=0.000 eV' 500 501 def distribution(self, eig_n, fermi_level): 502 f_n = np.zeros_like(eig_n) 503 f_n[eig_n < fermi_level] = 1.0 504 f_n[eig_n == fermi_level] = 0.5 505 return f_n, np.zeros_like(eig_n), np.zeros_like(eig_n) 506 507 def _calculate(self, 508 nelectrons, 509 eig_qn, 510 weight_q, 511 f_qn, 512 fermi_level_guess=nan): 513 eig_kn, weight_k, nkpts_r = collect_eigelvalues(eig_qn, weight_q, 514 self.bd, self.kpt_comm) 515 516 if eig_kn.size != 0: 517 # Try to use integer weights (avoid round-off errors): 518 N = int(round(1 / min(weight_k))) 519 w_k = (weight_k * N).round().astype(int) 520 if abs(w_k - N * weight_k).max() > 1e-10: 521 # Did not work. Use original fractional weights: 522 w_k = weight_k 523 N = 1 524 525 f_kn = np.zeros_like(eig_kn) 526 f_m = f_kn.ravel() 527 w_kn = np.empty_like(eig_kn, dtype=w_k.dtype) 528 w_kn[:] = w_k[:, np.newaxis] 529 eig_m = eig_kn.ravel() 530 w_m = w_kn.ravel() 531 m_i = eig_m.argsort() 532 w_i = w_m[m_i] 533 sum_i = np.add.accumulate(w_i) 534 filled_i = (sum_i <= nelectrons * N) 535 i = sum(filled_i) 536 f_m[m_i[:i]] = 1.0 537 if i == len(m_i): 538 fermi_level = inf 539 else: 540 extra = nelectrons * N - (sum_i[i - 1] if i > 0 else 0.0) 541 if extra > 0: 542 assert extra <= w_i[i] 543 f_m[m_i[i]] = extra / w_i[i] 544 fermi_level = eig_m[m_i[i]] 545 else: 546 fermi_level = (eig_m[m_i[i]] + eig_m[m_i[i - 1]]) / 2 547 else: 548 f_kn = None 549 fermi_level = nan 550 551 distribute_occupation_numbers(f_kn, f_qn, nkpts_r, 552 self.bd, self.kpt_comm) 553 554 if self.kpt_comm.rank == 0: 555 fermi_level = broadcast_float(fermi_level, self.bd.comm) 556 fermi_level = broadcast_float(fermi_level, self.kpt_comm) 557 558 e_entropy = 0.0 559 return fermi_level, e_entropy 560 561 562class FixedOccupationNumbers(OccupationNumberCalculator): 563 extrapolate_factor = 0.0 564 565 def __init__(self, numbers, parallel_layout: ParallelLayout = None): 566 """Fixed occupation numbers. 567 568 f_sn: ndarray, shape=(nspins, nbands) 569 Occupation numbers (in the range from 0 to 1) 570 571 Example (excited state with 4 electrons):: 572 573 occ = FixedOccupationNumbers([[1, 0, 1, 0], [1, 1, 0, 0]]) 574 575 """ 576 OccupationNumberCalculator.__init__(self, parallel_layout) 577 self.f_sn = np.array(numbers) 578 579 def _calculate(self, 580 nelectrons, 581 eig_qn, 582 weight_q, 583 f_qn, 584 fermi_level_guess=nan): 585 if self.bd.nbands == self.f_sn.shape[1]: 586 for q, f_n in enumerate(f_qn): 587 s = q % len(self.f_sn) 588 self.bd.distribute(self.f_sn[s], f_n) 589 else: 590 # Non-collinear calculation: 591 self.bd.distribute(self.f_sn.T.flatten().copy(), f_qn[0]) 592 593 return inf, 0.0 594 595 def todict(self): 596 return {'name': 'fixed', 'numbers': self.f_sn} 597 598 599def FixedOccupations(f_sn): 600 warnings.warn( 601 "Please use occupations={'name': 'fixed', 'numbers': ...} instead.") 602 if len(f_sn) == 1: 603 f_sn = np.array(f_sn) / 2 604 return {'name': 'fixed', 'numbers': f_sn} 605 606 607class ThomasFermiOccupations(OccupationNumberCalculator): 608 name = 'orbital-free' 609 extrapolate_factor = 0.0 610 611 def _calculate(self, 612 nelectrons, 613 eig_qn, 614 weight_q, 615 f_qn, 616 fermi_level_guess=nan): 617 assert len(f_qn) == 1 618 f_qn[0][:] = [nelectrons] 619 return inf, 0.0 620 621 622def create_occ_calc(dct: Dict[str, Any], 623 *, 624 parallel_layout: ParallelLayout = None, 625 fixed_magmom_value=None, 626 rcell=None, 627 monkhorst_pack_size=None, 628 bz2ibzmap=None, 629 ) -> OccupationNumberCalculator: 630 """Surprise: Create occupation-number object. 631 632 The unit of width is eV and name must be one of: 633 634 * 'fermi-dirac' 635 * 'marzari-vanderbilt' 636 * 'methfessel-paxton' 637 * 'fixed' 638 * 'tetrahedron-method' 639 * 'improved-tetrahedron-method' 640 * 'orbital-free' 641 642 >>> occ = create_occ_calc({'width': 0.0}) 643 >>> occ.calculate(nelectrons=3, 644 ... eigenvalues=[[0, 1, 2], [0, 2, 3]], 645 ... weights=[1, 1]) 646 (array([[1., 1., 0.], 647 [1., 0., 0.]]), [1.5], 0.0) 648 """ 649 kwargs = dct.copy() 650 fix_the_magnetic_moment = kwargs.pop('fixmagmom', False) 651 name = kwargs.pop('name', '') 652 kwargs['parallel_layout'] = parallel_layout 653 654 if name == 'unknown': 655 return OccupationNumberCalculator(**kwargs) 656 657 occ: OccupationNumberCalculator 658 659 if kwargs.get('width') == 0.0: 660 del kwargs['width'] 661 occ = ZeroWidth(**kwargs) 662 elif name == 'methfessel-paxton': 663 occ = MethfesselPaxtonCalculator(**kwargs) 664 elif name == 'fermi-dirac': 665 occ = FermiDiracCalculator(**kwargs) 666 elif name == 'marzari-vanderbilt': 667 occ = MarzariVanderbiltCalculator(**kwargs) 668 elif name in {'tetrahedron-method', 'improved-tetrahedron-method'}: 669 from gpaw.tetrahedron import TetrahedronMethod 670 occ = TetrahedronMethod(rcell, 671 monkhorst_pack_size, 672 name == 'improved-tetrahedron-method', 673 bz2ibzmap, 674 **kwargs) 675 elif name == 'orbital-free': 676 return ThomasFermiOccupations(**kwargs) 677 elif name == 'fixed': 678 return FixedOccupationNumbers(**kwargs) 679 else: 680 raise ValueError(f'Unknown occupation number object name: {name}') 681 682 if fix_the_magnetic_moment: 683 occ = FixMagneticMomentOccupationNumberCalculator( 684 occ, fixed_magmom_value) 685 686 return occ 687 688 689def occupation_numbers(occ, eig_skn, weight_k, nelectrons): 690 """Calculate occupation numbers from eigenvalues in eV (**deprecated**). 691 692 occ: dict 693 Example: {'name': 'fermi-dirac', 'width': 0.05} (width in eV). 694 eps_skn: ndarray, shape=(nspins, nibzkpts, nbands) 695 Eigenvalues. 696 weight_k: ndarray, shape=(nibzkpts,) 697 Weights of k-points in IBZ (must sum to 1). 698 nelectrons: int or float 699 Number of electrons. 700 701 Returns a tuple containing: 702 703 * f_skn (sums to nelectrons) 704 * fermi-level [Hartree] 705 * magnetic moment 706 * entropy as -S*T [Hartree] 707 """ 708 709 warnings.warn('Please use one of the OccupationNumbers implementations', 710 DeprecationWarning) 711 occ = create_occ_calc(occ) 712 f_kn, (fermi_level,), e_entropy = occ.calculate( 713 nelectrons * len(eig_skn) / 2, 714 [eig_n for eig_kn in eig_skn for eig_n in eig_kn], 715 list(weight_k) * len(eig_skn)) 716 717 f_kn *= np.array(weight_k)[:, np.newaxis] 718 719 if len(eig_skn) == 1: 720 f_skn = np.array([f_kn]) * 2 721 e_entropy *= 2 722 magmom = 0.0 723 else: 724 f_skn = np.array([f_kn[::2], f_kn[1::2]]) 725 f1, f2 = f_skn.sum(axis=(1, 2)) 726 magmom = f1 - f2 727 728 return f_skn, fermi_level * Ha, magmom, e_entropy * Ha 729 730 731def FermiDirac(width, fixmagmom=False): 732 return dict(name='fermi-dirac', width=width, fixmagmom=fixmagmom) 733 734 735def MarzariVanderbilt(width, fixmagmom=False): 736 return dict(name='marzari-vanderbilt', width=width, fixmagmom=fixmagmom) 737 738 739def MethfesselPaxton(width, order=0, fixmagmom=False): 740 return dict(name='methfessel-paxton', width=width, order=order, 741 fixmagmom=fixmagmom) 742