1# coding: utf-8 2 3""" 4This module provides classes to run and analyze boltztrap on pymatgen band 5structure objects. Boltztrap is a software interpolating band structures and 6computing materials properties from this band structure using Boltzmann 7semi-classical transport theory. 8 9Boltztrap has been developed by Georg Madsen. 10 11http://www.icams.de/content/research/software-development/boltztrap/ 12 13You need version 1.2.3 or higher 14 15References are:: 16 17 Madsen, G. K. H., and Singh, D. J. (2006). 18 BoltzTraP. A code for calculating band-structure dependent quantities. 19 Computer Physics Communications, 175, 67-71 20""" 21 22import logging 23import math 24import os 25import subprocess 26import tempfile 27import time 28 29import numpy as np 30from monty.dev import requires 31from monty.json import MSONable, jsanitize 32from monty.os import cd 33from monty.os.path import which 34from scipy import constants 35from scipy.spatial import distance 36 37from pymatgen.core.lattice import Lattice 38from pymatgen.core.units import Energy, Length 39from pymatgen.electronic_structure.bandstructure import BandStructureSymmLine, Kpoint 40from pymatgen.electronic_structure.core import Orbital 41from pymatgen.electronic_structure.dos import CompleteDos, Dos, Spin 42from pymatgen.symmetry.analyzer import SpacegroupAnalyzer 43from pymatgen.symmetry.bandstructure import HighSymmKpath 44 45__author__ = "Geoffroy Hautier, Zachary Gibbs, Francesco Ricci, Anubhav Jain" 46__copyright__ = "Copyright 2013, The Materials Project" 47__version__ = "1.1" 48__maintainer__ = "Geoffroy Hautier" 49__email__ = "geoffroy@uclouvain.be" 50__status__ = "Development" 51__date__ = "August 23, 2013" 52 53 54class BoltztrapRunner(MSONable): 55 """ 56 This class is used to run Boltztrap on a band structure object. 57 """ 58 59 @requires( 60 which("x_trans"), 61 "BoltztrapRunner requires the executables 'x_trans' to be in " 62 "the path. Please download the Boltztrap at http://" 63 "www.icams.de/content/research/software-development/boltztrap/ " 64 "and follow the instructions in the README to compile " 65 "Bolztrap accordingly. Then add x_trans to your path", 66 ) 67 def __init__( 68 self, 69 bs, 70 nelec, 71 dos_type="HISTO", 72 energy_grid=0.005, 73 lpfac=10, 74 run_type="BOLTZ", 75 band_nb=None, 76 tauref=0, 77 tauexp=0, 78 tauen=0, 79 soc=False, 80 doping=None, 81 energy_span_around_fermi=1.5, 82 scissor=0.0, 83 kpt_line=None, 84 spin=None, 85 cond_band=False, 86 tmax=1300, 87 tgrid=50, 88 symprec=1e-3, 89 cb_cut=10, 90 timeout=7200, 91 ): 92 """ 93 Args: 94 bs: 95 A band structure object 96 nelec: 97 the number of electrons 98 dos_type: 99 two options for the band structure integration: "HISTO" 100 (histogram) or "TETRA" using the tetrahedon method. TETRA 101 typically gives better results (especially for DOSes) 102 but takes more time 103 energy_grid: 104 the energy steps used for the integration (eV) 105 lpfac: 106 the number of interpolation points in the real space. By 107 default 10 gives 10 time more points in the real space than 108 the number of kpoints given in reciprocal space 109 run_type: 110 type of boltztrap usage. by default 111 - BOLTZ: (default) compute transport coefficients 112 - BANDS: interpolate all bands contained in the energy range 113 specified in energy_span_around_fermi variable, along specified 114 k-points 115 - DOS: compute total and partial dos (custom BoltzTraP code 116 needed!) 117 - FERMI: compute fermi surface or more correctly to 118 get certain bands interpolated 119 band_nb: 120 indicates a band number. Used for Fermi Surface interpolation 121 (run_type="FERMI") 122 spin: 123 specific spin component (1: up, -1: down) of the band selected 124 in FERMI mode (mandatory). 125 cond_band: 126 if a conduction band is specified in FERMI mode, 127 set this variable as True 128 tauref: 129 reference relaxation time. Only set to a value different than 130 zero if we want to model beyond the constant relaxation time. 131 tauexp: 132 exponent for the energy in the non-constant relaxation time 133 approach 134 tauen: 135 reference energy for the non-constant relaxation time approach 136 soc: 137 results from spin-orbit coupling (soc) computations give 138 typically non-polarized (no spin up or down) results but single 139 electron occupations. If the band structure comes from a soc 140 computation, you should set soc to True (default False) 141 doping: 142 the fixed doping levels you want to compute. Boltztrap provides 143 both transport values depending on electron chemical potential 144 (fermi energy) and for a series of fixed carrier 145 concentrations. By default, this is set to 1e16 to 1e22 in 146 increments of factors of 10. 147 energy_span_around_fermi: 148 usually the interpolation is not needed on the entire energy 149 range but on a specific range around the fermi level. 150 This energy gives this range in eV. by default it is 1.5 eV. 151 If DOS or BANDS type are selected, this range is automatically 152 set to cover the entire energy range. 153 scissor: 154 scissor to apply to the band gap (eV). This applies a scissor 155 operation moving the band edges without changing the band 156 shape. This is useful to correct the often underestimated band 157 gap in DFT. Default is 0.0 (no scissor) 158 kpt_line: 159 list of fractional coordinates of kpoints as arrays or list of 160 Kpoint objects for BANDS mode calculation (standard path of 161 high symmetry k-points is automatically set as default) 162 tmax: 163 Maximum temperature (K) for calculation (default=1300) 164 tgrid: 165 Temperature interval for calculation (default=50) 166 symprec: 1e-3 is the default in pymatgen. If the kmesh has been 167 generated using a different symprec, it has to be specified 168 to avoid a "factorization error" in BoltzTraP calculation. 169 If a kmesh that spans the whole Brillouin zone has been used, 170 or to disable all the symmetries, set symprec to None. 171 cb_cut: by default 10% of the highest conduction bands are 172 removed because they are often not accurate. 173 Tune cb_cut to change the percentage (0-100) of bands 174 that are removed. 175 timeout: overall time limit (in seconds): mainly to avoid infinite 176 loop when trying to find Fermi levels. 177 """ 178 self.lpfac = lpfac 179 self._bs = bs 180 self._nelec = nelec 181 self.dos_type = dos_type 182 self.energy_grid = energy_grid 183 self.error = [] 184 self.run_type = run_type 185 self.band_nb = band_nb 186 self.spin = spin 187 self.cond_band = cond_band 188 self.tauref = tauref 189 self.tauexp = tauexp 190 self.tauen = tauen 191 self.soc = soc 192 self.kpt_line = kpt_line 193 self.cb_cut = cb_cut / 100.0 194 if isinstance(doping, list) and len(doping) > 0: 195 self.doping = doping 196 else: 197 self.doping = [] 198 for d in [1e16, 1e17, 1e18, 1e19, 1e20, 1e21]: 199 self.doping.extend([1 * d, 2.5 * d, 5 * d, 7.5 * d]) 200 self.doping.append(1e22) 201 self.energy_span_around_fermi = energy_span_around_fermi 202 self.scissor = scissor 203 self.tmax = tmax 204 self.tgrid = tgrid 205 self._symprec = symprec 206 if self.run_type in ("DOS", "BANDS"): 207 self._auto_set_energy_range() 208 self.timeout = timeout 209 self.start_time = time.time() 210 211 def _auto_set_energy_range(self): 212 """ 213 automatically determine the energy range as min/max eigenvalue 214 minus/plus the buffer_in_ev 215 """ 216 emins = [min([e_k[0] for e_k in self._bs.bands[Spin.up]])] 217 emaxs = [max([e_k[0] for e_k in self._bs.bands[Spin.up]])] 218 219 if self._bs.is_spin_polarized: 220 emins.append(min([e_k[0] for e_k in self._bs.bands[Spin.down]])) 221 222 emaxs.append(max([e_k[0] for e_k in self._bs.bands[Spin.down]])) 223 224 min_eigenval = Energy(min(emins) - self._bs.efermi, "eV").to("Ry") 225 max_eigenval = Energy(max(emaxs) - self._bs.efermi, "eV").to("Ry") 226 227 # set energy range to buffer around min/max EV 228 # buffer does not increase CPU time but will help get equal 229 # energies for spin up/down for band structure 230 const = Energy(2, "eV").to("Ry") 231 self._ll = min_eigenval - const 232 self._hl = max_eigenval + const 233 234 en_range = Energy(max((abs(self._ll), abs(self._hl))), "Ry").to("eV") 235 236 self.energy_span_around_fermi = en_range * 1.01 237 print("energy_span_around_fermi = ", self.energy_span_around_fermi) 238 239 @property 240 def bs(self): 241 """ 242 :return: The BandStructure 243 """ 244 return self._bs 245 246 @property 247 def nelec(self): 248 """ 249 :return: Number of electrons 250 """ 251 return self._nelec 252 253 def write_energy(self, output_file): 254 """ 255 Writes the energy to an output file. 256 257 :param output_file: Filename 258 """ 259 with open(output_file, "w") as f: 260 f.write("test\n") 261 f.write("{}\n".format(len(self._bs.kpoints))) 262 263 if self.run_type == "FERMI": 264 sign = -1.0 if self.cond_band else 1.0 265 for i, kpt in enumerate(self._bs.kpoints): 266 eigs = [] 267 eigs.append( 268 Energy( 269 self._bs.bands[Spin(self.spin)][self.band_nb][i] - self._bs.efermi, 270 "eV", 271 ).to("Ry") 272 ) 273 f.write( 274 "%12.8f %12.8f %12.8f %d\n" 275 % ( 276 kpt.frac_coords[0], 277 kpt.frac_coords[1], 278 kpt.frac_coords[2], 279 len(eigs), 280 ) 281 ) 282 for e in eigs: 283 f.write("%18.8f\n" % (sign * float(e))) 284 285 else: 286 for i, kpt in enumerate(self._bs.kpoints): 287 eigs = [] 288 if self.run_type == "DOS": 289 spin_lst = [self.spin] 290 else: 291 spin_lst = self._bs.bands 292 293 for spin in spin_lst: 294 # use 90% of bottom bands since highest eigenvalues 295 # are usually incorrect 296 # ask Geoffroy Hautier for more details 297 nb_bands = int(math.floor(self._bs.nb_bands * (1 - self.cb_cut))) 298 for j in range(nb_bands): 299 eigs.append( 300 Energy( 301 self._bs.bands[Spin(spin)][j][i] - self._bs.efermi, 302 "eV", 303 ).to("Ry") 304 ) 305 eigs.sort() 306 307 if self.run_type == "DOS" and self._bs.is_spin_polarized: 308 eigs.insert(0, self._ll) 309 eigs.append(self._hl) 310 311 f.write( 312 "%12.8f %12.8f %12.8f %d\n" 313 % ( 314 kpt.frac_coords[0], 315 kpt.frac_coords[1], 316 kpt.frac_coords[2], 317 len(eigs), 318 ) 319 ) 320 321 for e in eigs: 322 f.write("%18.8f\n" % (float(e))) 323 324 def write_struct(self, output_file): 325 """ 326 Writes the structure to an output file. 327 328 :param output_file: Filename 329 """ 330 if self._symprec is not None: 331 sym = SpacegroupAnalyzer(self._bs.structure, symprec=self._symprec) 332 elif self._symprec is None: 333 pass 334 335 with open(output_file, "w") as f: 336 if self._symprec is not None: 337 f.write( 338 "{} {}\n".format( 339 self._bs.structure.composition.formula, 340 sym.get_space_group_symbol(), 341 ) 342 ) 343 elif self._symprec is None: 344 f.write("{} {}\n".format(self._bs.structure.composition.formula, "symmetries disabled")) 345 346 f.write( 347 "{}\n".format( 348 "\n".join( 349 [ 350 " ".join(["%.5f" % Length(i, "ang").to("bohr") for i in row]) 351 for row in self._bs.structure.lattice.matrix 352 ] 353 ) 354 ) 355 ) 356 357 if self._symprec is not None: 358 ops = sym.get_symmetry_dataset()["rotations"] 359 elif self._symprec is None: 360 ops = [[[1, 0, 0], [0, 1, 0], [0, 0, 1]]] 361 f.write("{}\n".format(len(ops))) 362 363 for c in ops: 364 for row in c: 365 f.write("{}\n".format(" ".join(str(i) for i in row))) 366 367 def write_def(self, output_file): 368 """ 369 Writes the def to an output file. 370 371 :param output_file: Filename 372 """ 373 # This function is useless in std version of BoltzTraP code 374 # because x_trans script overwrite BoltzTraP.def 375 with open(output_file, "w") as f: 376 so = "" 377 if self._bs.is_spin_polarized or self.soc: 378 so = "so" 379 f.write( 380 "5, 'boltztrap.intrans', 'old', 'formatted',0\n" 381 + "6,'boltztrap.outputtrans', 'unknown', " 382 "'formatted',0\n" 383 + "20,'boltztrap.struct', 'old', 'formatted',0\n" 384 + "10,'boltztrap.energy" 385 + so 386 + "', 'old', " 387 "'formatted',0\n" + "48,'boltztrap.engre', 'unknown', " 388 "'unformatted',0\n" + "49,'boltztrap.transdos', 'unknown', " 389 "'formatted',0\n" + "50,'boltztrap.sigxx', 'unknown', 'formatted'," 390 "0\n" + "51,'boltztrap.sigxxx', 'unknown', 'formatted'," 391 "0\n" + "21,'boltztrap.trace', 'unknown', " 392 "'formatted',0\n" + "22,'boltztrap.condtens', 'unknown', " 393 "'formatted',0\n" + "24,'boltztrap.halltens', 'unknown', " 394 "'formatted',0\n" + "30,'boltztrap_BZ.cube', 'unknown', " 395 "'formatted',0\n" 396 ) 397 398 def write_proj(self, output_file_proj, output_file_def): 399 """ 400 Writes the projections to an output file. 401 402 :param output_file: Filename 403 """ 404 # This function is useless in std version of BoltzTraP code 405 # because x_trans script overwrite BoltzTraP.def 406 for oi, o in enumerate(Orbital): 407 for site_nb in range(0, len(self._bs.structure.sites)): 408 if oi < len(self._bs.projections[Spin.up][0][0]): 409 with open(output_file_proj + "_" + str(site_nb) + "_" + str(o), "w") as f: 410 f.write(self._bs.structure.composition.formula + "\n") 411 f.write(str(len(self._bs.kpoints)) + "\n") 412 for i, kpt in enumerate(self._bs.kpoints): 413 tmp_proj = [] 414 for j in range(int(math.floor(self._bs.nb_bands * (1 - self.cb_cut)))): 415 tmp_proj.append(self._bs.projections[Spin(self.spin)][j][i][oi][site_nb]) 416 # TODO deal with the sorting going on at 417 # the energy level!!! 418 # tmp_proj.sort() 419 420 if self.run_type == "DOS" and self._bs.is_spin_polarized: 421 tmp_proj.insert(0, self._ll) 422 tmp_proj.append(self._hl) 423 424 f.write( 425 "%12.8f %12.8f %12.8f %d\n" 426 % ( 427 kpt.frac_coords[0], 428 kpt.frac_coords[1], 429 kpt.frac_coords[2], 430 len(tmp_proj), 431 ) 432 ) 433 for t in tmp_proj: 434 f.write("%18.8f\n" % float(t)) 435 with open(output_file_def, "w") as f: 436 so = "" 437 if self._bs.is_spin_polarized: 438 so = "so" 439 f.write( 440 "5, 'boltztrap.intrans', 'old', 'formatted',0\n" 441 + "6,'boltztrap.outputtrans', 'unknown', " 442 "'formatted',0\n" 443 + "20,'boltztrap.struct', 'old', 'formatted',0\n" 444 + "10,'boltztrap.energy" 445 + so 446 + "', 'old', " 447 "'formatted',0\n" + "48,'boltztrap.engre', 'unknown', " 448 "'unformatted',0\n" + "49,'boltztrap.transdos', 'unknown', " 449 "'formatted',0\n" + "50,'boltztrap.sigxx', 'unknown', 'formatted'," 450 "0\n" + "51,'boltztrap.sigxxx', 'unknown', 'formatted'," 451 "0\n" + "21,'boltztrap.trace', 'unknown', " 452 "'formatted',0\n" + "22,'boltztrap.condtens', 'unknown', " 453 "'formatted',0\n" + "24,'boltztrap.halltens', 'unknown', " 454 "'formatted',0\n" + "30,'boltztrap_BZ.cube', 'unknown', " 455 "'formatted',0\n" 456 ) 457 i = 1000 458 for oi, o in enumerate(Orbital): 459 for site_nb in range(0, len(self._bs.structure.sites)): 460 if oi < len(self._bs.projections[Spin.up][0][0]): 461 f.write( 462 str(i) 463 + ",'" 464 + "boltztrap.proj_" 465 + str(site_nb) 466 + "_" 467 + str(o.name) 468 + "' 'old', 'formatted',0\n" 469 ) 470 i += 1 471 472 def write_intrans(self, output_file): 473 """ 474 Writes the intrans to an output file. 475 476 :param output_file: Filename 477 """ 478 setgap = 1 if self.scissor > 0.0001 else 0 479 480 if self.run_type in ("BOLTZ", "DOS"): 481 with open(output_file, "w") as fout: 482 fout.write("GENE # use generic interface\n") 483 fout.write( 484 "1 0 %d %f # iskip (not presently used) idebug " 485 "setgap shiftgap \n" % (setgap, Energy(self.scissor, "eV").to("Ry")) 486 ) 487 fout.write( 488 "0.0 %f %f %6.1f # Fermilevel (Ry),energygrid,energy " 489 "span around Fermilevel, number of electrons\n" 490 % ( 491 Energy(self.energy_grid, "eV").to("Ry"), 492 Energy(self.energy_span_around_fermi, "eV").to("Ry"), 493 self._nelec, 494 ) 495 ) 496 fout.write("CALC # CALC (calculate expansion " "coeff), NOCALC read from file\n") 497 fout.write("%d # lpfac, number of latt-points " "per k-point\n" % self.lpfac) 498 fout.write("%s # run mode (only BOLTZ is " "supported)\n" % self.run_type) 499 fout.write(".15 # (efcut) energy range of " "chemical potential\n") 500 fout.write("{} {} # Tmax, temperature grid\n".format(self.tmax, self.tgrid)) 501 fout.write("-1. # energyrange of bands given DOS output sig_xxx and " "dos_xxx (xxx is band number)\n") 502 fout.write(self.dos_type + "\n") # e.g., HISTO or TETRA 503 fout.write("{} {} {} 0 0 0\n".format(self.tauref, self.tauexp, self.tauen)) 504 fout.write("{}\n".format(2 * len(self.doping))) 505 506 for d in self.doping: 507 fout.write(str(d) + "\n") 508 for d in self.doping: 509 fout.write(str(-d) + "\n") 510 511 elif self.run_type == "FERMI": 512 with open(output_file, "w") as fout: 513 fout.write("GENE # use generic interface\n") 514 fout.write("1 0 0 0.0 # iskip (not presently used) idebug " "setgap shiftgap \n") 515 fout.write( 516 "0.0 %f 0.1 %6.1f # Fermilevel (Ry),energygrid," 517 "energy span around Fermilevel, " 518 "number of electrons\n" % (Energy(self.energy_grid, "eV").to("Ry"), self._nelec) 519 ) 520 fout.write("CALC # CALC (calculate expansion " "coeff), NOCALC read from file\n") 521 fout.write("%d # lpfac, number of latt-points " "per k-point\n" % self.lpfac) 522 fout.write("FERMI # run mode (only BOLTZ is " "supported)\n") 523 fout.write( 524 str(1) 525 + " # actual band selected: " 526 + str(self.band_nb + 1) 527 + " spin: " 528 + str(self.spin) 529 ) 530 531 elif self.run_type == "BANDS": 532 if self.kpt_line is None: 533 kpath = HighSymmKpath(self._bs.structure) 534 self.kpt_line = [ 535 Kpoint(k, self._bs.structure.lattice) for k in kpath.get_kpoints(coords_are_cartesian=False)[0] 536 ] 537 self.kpt_line = [kp.frac_coords for kp in self.kpt_line] 538 elif isinstance(self.kpt_line[0], Kpoint): 539 self.kpt_line = [kp.frac_coords for kp in self.kpt_line] 540 541 with open(output_file, "w") as fout: 542 fout.write("GENE # use generic interface\n") 543 fout.write( 544 "1 0 %d %f # iskip (not presently used) idebug " 545 "setgap shiftgap \n" % (setgap, Energy(self.scissor, "eV").to("Ry")) 546 ) 547 fout.write( 548 "0.0 %f %f %6.1f # Fermilevel (Ry),energygrid,energy " 549 "span around Fermilevel, " 550 "number of electrons\n" 551 % ( 552 Energy(self.energy_grid, "eV").to("Ry"), 553 Energy(self.energy_span_around_fermi, "eV").to("Ry"), 554 self._nelec, 555 ) 556 ) 557 fout.write("CALC # CALC (calculate expansion " "coeff), NOCALC read from file\n") 558 fout.write("%d # lpfac, number of latt-points " "per k-point\n" % self.lpfac) 559 fout.write("BANDS # run mode (only BOLTZ is " "supported)\n") 560 fout.write("P " + str(len(self.kpt_line)) + "\n") 561 for kp in self.kpt_line: 562 fout.writelines([str(k) + " " for k in kp]) 563 fout.write("\n") 564 565 def write_input(self, output_dir): 566 """ 567 Writes the input files. 568 569 :param output_dir: Directory to write the input files. 570 """ 571 if self._bs.is_spin_polarized or self.soc: 572 self.write_energy(os.path.join(output_dir, "boltztrap.energyso")) 573 else: 574 self.write_energy(os.path.join(output_dir, "boltztrap.energy")) 575 576 self.write_struct(os.path.join(output_dir, "boltztrap.struct")) 577 self.write_intrans(os.path.join(output_dir, "boltztrap.intrans")) 578 self.write_def(os.path.join(output_dir, "BoltzTraP.def")) 579 580 if len(self.bs.projections) != 0 and self.run_type == "DOS": 581 self.write_proj( 582 os.path.join(output_dir, "boltztrap.proj"), 583 os.path.join(output_dir, "BoltzTraP.def"), 584 ) 585 586 def run( 587 self, 588 path_dir=None, 589 convergence=True, 590 write_input=True, 591 clear_dir=False, 592 max_lpfac=150, 593 min_egrid=0.00005, 594 ): 595 """ 596 Write inputs (optional), run BoltzTraP, and ensure 597 convergence (optional) 598 Args: 599 path_dir (str): directory in which to run BoltzTraP 600 convergence (bool): whether to check convergence and make 601 corrections if needed 602 write_input: (bool) whether to write input files before the run 603 (required for convergence mode) 604 clear_dir: (bool) whether to remove all files in the path_dir 605 before starting 606 max_lpfac: (float) maximum lpfac value to try before reducing egrid 607 in convergence mode 608 min_egrid: (float) minimum egrid value to try before giving up in 609 convergence mode 610 611 Returns: 612 613 """ 614 615 # TODO: consider making this a part of custodian rather than pymatgen 616 # A lot of this functionality (scratch dirs, handlers, monitors) 617 # is built into custodian framework 618 619 if convergence and not write_input: 620 raise ValueError("Convergence mode requires write_input to be " "true") 621 622 if self.run_type in ("BANDS", "DOS", "FERMI"): 623 convergence = False 624 if self.lpfac > max_lpfac: 625 max_lpfac = self.lpfac 626 627 if self.run_type == "BANDS" and self.bs.is_spin_polarized: 628 print( 629 "Reminder: for run_type " + str(self.run_type) + ", spin component are not separated! " 630 "(you have a spin polarized band structure)" 631 ) 632 633 if self.run_type in ("FERMI", "DOS") and self.spin is None: 634 if self.bs.is_spin_polarized: 635 raise BoltztrapError("Spin parameter must be specified for spin polarized " "band structures!") 636 self.spin = 1 637 638 dir_bz_name = "boltztrap" 639 if path_dir is None: 640 temp_dir = tempfile.mkdtemp() 641 path_dir = os.path.join(temp_dir, dir_bz_name) 642 else: 643 path_dir = os.path.abspath(os.path.join(path_dir, dir_bz_name)) 644 645 if not os.path.exists(path_dir): 646 os.mkdir(path_dir) 647 elif clear_dir: 648 for c in os.listdir(path_dir): 649 os.remove(os.path.join(path_dir, c)) 650 651 FORMAT = "%(message)s" 652 logging.basicConfig( 653 level=logging.INFO, 654 format=FORMAT, 655 filename=os.path.join(path_dir, "../boltztrap.out"), 656 ) 657 658 with cd(path_dir): 659 lpfac_start = self.lpfac 660 converged = False 661 662 while self.energy_grid >= min_egrid and not converged: 663 self.lpfac = lpfac_start 664 if time.time() - self.start_time > self.timeout: 665 raise BoltztrapError("no doping convergence after timeout " "of {} s".format(self.timeout)) 666 667 logging.info("lpfac, energy_grid: {} {}".format(self.lpfac, self.energy_grid)) 668 669 while self.lpfac <= max_lpfac and not converged: 670 if time.time() - self.start_time > self.timeout: 671 raise BoltztrapError("no doping convergence after " "timeout of {} s".format(self.timeout)) 672 673 if write_input: 674 self.write_input(path_dir) 675 676 bt_exe = ["x_trans", "BoltzTraP"] 677 if self._bs.is_spin_polarized or self.soc: 678 bt_exe.append("-so") 679 680 with subprocess.Popen( 681 bt_exe, 682 stdout=subprocess.PIPE, 683 stdin=subprocess.PIPE, 684 stderr=subprocess.PIPE, 685 ) as p: 686 p.wait() 687 688 for c in p.communicate(): 689 logging.info(c.decode()) 690 if "error in factorization" in c.decode(): 691 raise BoltztrapError("error in factorization") 692 693 warning = "" 694 695 with open(os.path.join(path_dir, dir_bz_name + ".outputtrans")) as f: 696 for l in f: 697 if "Option unknown" in l: 698 raise BoltztrapError("DOS mode needs a custom version of " "BoltzTraP code is needed") 699 if "WARNING" in l: 700 warning = l 701 break 702 if "Error - Fermi level was not found" in l: 703 warning = l 704 break 705 706 if not warning and convergence: 707 # check convergence for warning 708 analyzer = BoltztrapAnalyzer.from_files(path_dir) 709 for doping in ["n", "p"]: 710 for c in analyzer.mu_doping[doping]: 711 if len(analyzer.mu_doping[doping][c]) != len(analyzer.doping[doping]): 712 warning = "length of mu_doping array is " "incorrect" 713 break 714 715 if ( 716 doping == "p" 717 and sorted(analyzer.mu_doping[doping][c], reverse=True) 718 != analyzer.mu_doping[doping][c] 719 ): 720 warning = "sorting of mu_doping array " "incorrect for p-type" 721 break 722 723 # ensure n-type doping sorted correctly 724 if ( 725 doping == "n" 726 and sorted(analyzer.mu_doping[doping][c]) != analyzer.mu_doping[doping][c] 727 ): 728 warning = "sorting of mu_doping array " "incorrect for n-type" 729 break 730 731 if warning: 732 self.lpfac += 10 733 logging.warn("Warning detected: {}! Increase lpfac to " "{}".format(warning, self.lpfac)) 734 735 else: 736 converged = True 737 738 if not converged: 739 self.energy_grid /= 10 740 logging.info("Could not converge with max lpfac; " "Decrease egrid to {}".format(self.energy_grid)) 741 742 if not converged: 743 raise BoltztrapError( 744 "Doping convergence not reached with lpfac=" 745 + str(self.lpfac) 746 + ", energy_grid=" 747 + str(self.energy_grid) 748 ) 749 750 return path_dir 751 752 def as_dict(self): 753 """ 754 :return: MSONable dict 755 """ 756 results = { 757 "@module": self.__class__.__module__, 758 "@class": self.__class__.__name__, 759 "lpfac": self.lpfac, 760 "bs": self.bs.as_dict(), 761 "nelec": self._nelec, 762 "dos_type": self.dos_type, 763 "run_type": self.run_type, 764 "band_nb": self.band_nb, 765 "spin": self.spin, 766 "cond_band": self.cond_band, 767 "tauref": self.tauref, 768 "tauexp": self.tauexp, 769 "tauen": self.tauen, 770 "soc": self.soc, 771 "kpt_line": self.kpt_line, 772 "doping": self.doping, 773 "energy_span_around_fermi": self.energy_span_around_fermi, 774 "scissor": self.scissor, 775 "tmax": self.tmax, 776 "tgrid": self.tgrid, 777 "symprec": self._symprec, 778 } 779 return jsanitize(results) 780 781 782class BoltztrapError(Exception): 783 """ 784 Exception class for boltztrap. 785 Raised when the boltztrap gives an error 786 """ 787 788 pass 789 790 791class BoltztrapAnalyzer: 792 """ 793 Class used to store all the data from a boltztrap run 794 """ 795 796 def __init__( 797 self, 798 gap=None, 799 mu_steps=None, 800 cond=None, 801 seebeck=None, 802 kappa=None, 803 hall=None, 804 doping=None, 805 mu_doping=None, 806 seebeck_doping=None, 807 cond_doping=None, 808 kappa_doping=None, 809 hall_doping=None, 810 intrans=None, 811 dos=None, 812 dos_partial=None, 813 carrier_conc=None, 814 vol=None, 815 warning=None, 816 bz_bands=None, 817 bz_kpoints=None, 818 fermi_surface_data=None, 819 ): 820 """ 821 Constructor taking directly all the data generated by Boltztrap. You 822 won't probably use it directly but instead use the from_files and 823 from_dict methods. 824 825 Args: 826 gap: The gap after interpolation in eV 827 mu_steps: The steps of electron chemical potential (or Fermi 828 level) in eV. 829 cond: The electronic conductivity tensor divided by a constant 830 relaxation time (sigma/tau) at different temperature and 831 fermi levels. 832 The format is {temperature: [array of 3x3 tensors at each 833 fermi level in mu_steps]}. The units are 1/(Ohm*m*s). 834 seebeck: The Seebeck tensor at different temperatures and fermi 835 levels. The format is {temperature: [array of 3x3 tensors at 836 each fermi level in mu_steps]}. The units are V/K 837 kappa: The electronic thermal conductivity tensor divided by a 838 constant relaxation time (kappa/tau) at different temperature 839 and fermi levels. The format is {temperature: [array of 3x3 840 tensors at each fermi level in mu_steps]} 841 The units are W/(m*K*s) 842 hall: The hall tensor at different temperature and fermi levels 843 The format is {temperature: [array of 27 coefficients list at 844 each fermi level in mu_steps]} 845 The units are m^3/C 846 doping: The different doping levels that have been given to 847 Boltztrap. The format is {'p':[],'n':[]} with an array of 848 doping levels. The units are cm^-3 849 mu_doping: Gives the electron chemical potential (or Fermi level) 850 for a given set of doping. 851 Format is {'p':{temperature: [fermi levels],'n':{temperature: 852 [fermi levels]}} 853 the fermi level array is ordered according to the doping 854 levels in doping units for doping are in cm^-3 and for Fermi 855 level in eV 856 seebeck_doping: The Seebeck tensor at different temperatures and 857 doping levels. The format is {'p': {temperature: [Seebeck 858 tensors]}, 'n':{temperature: [Seebeck tensors]}} 859 The [Seebeck tensors] array is ordered according to the 860 doping levels in doping units for doping are in cm^-3 and for 861 Seebeck in V/K 862 cond_doping: The electronic conductivity tensor divided by a 863 constant relaxation time (sigma/tau) at different 864 temperatures and doping levels 865 The format is {'p':{temperature: [conductivity tensors]}, 866 'n':{temperature: [conductivity tensors]}} 867 The [conductivity tensors] array is ordered according to the 868 doping levels in doping units for doping are in cm^-3 and for 869 conductivity in 1/(Ohm*m*s) 870 kappa_doping: The thermal conductivity tensor divided by a constant 871 relaxation time (kappa/tau) at different temperatures and 872 doping levels. 873 The format is {'p':{temperature: [thermal conductivity 874 tensors]},'n':{temperature: [thermal conductivity tensors]}} 875 The [thermal conductivity tensors] array is ordered according 876 to the doping levels in doping units for doping are in cm^-3 877 and for thermal conductivity in W/(m*K*s) 878 hall_doping: The Hall tensor at different temperatures and doping 879 levels. 880 The format is {'p':{temperature: [Hall tensors]}, 881 'n':{temperature: [Hall tensors]}} 882 The [Hall tensors] array is ordered according to the doping 883 levels in doping and each Hall tensor is represented by a 27 884 coefficients list. 885 The units are m^3/C 886 intrans: a dictionary of inputs e.g. {"scissor": 0.0} 887 carrier_conc: The concentration of carriers in electron (or hole) 888 per unit cell 889 dos: The dos computed by Boltztrap given as a pymatgen Dos object 890 dos_partial: Data for the partial DOS projected on sites and 891 orbitals 892 vol: Volume of the unit cell in angstrom cube (A^3) 893 warning: string if BoltzTraP outputted a warning, else None 894 bz_bands: Data for interpolated bands on a k-point line 895 (run_type=BANDS) 896 bz_kpoints: k-point in reciprocal coordinates for interpolated 897 bands (run_type=BANDS) 898 fermi_surface_data: energy values in a 3D grid imported from the 899 output .cube file. 900 """ 901 self.gap = gap 902 self.mu_steps = mu_steps 903 self._cond = cond 904 self._seebeck = seebeck 905 self._kappa = kappa 906 self._hall = hall 907 self.warning = warning 908 self.doping = doping 909 self.mu_doping = mu_doping 910 self._seebeck_doping = seebeck_doping 911 self._cond_doping = cond_doping 912 self._kappa_doping = kappa_doping 913 self._hall_doping = hall_doping 914 self.intrans = intrans 915 self._carrier_conc = carrier_conc 916 self.dos = dos 917 self.vol = vol 918 self._dos_partial = dos_partial 919 self._bz_bands = bz_bands 920 self._bz_kpoints = bz_kpoints 921 self.fermi_surface_data = fermi_surface_data 922 923 def get_symm_bands(self, structure, efermi, kpt_line=None, labels_dict=None): 924 """ 925 Function useful to read bands from Boltztrap output and get a 926 BandStructureSymmLine object comparable with that one from a DFT 927 calculation (if the same kpt_line is provided). Default kpt_line 928 and labels_dict is the standard path of high symmetry k-point for 929 the specified structure. They could be extracted from the 930 BandStructureSymmLine object that you want to compare with. efermi 931 variable must be specified to create the BandStructureSymmLine 932 object (usually it comes from DFT or Boltztrap calc) 933 """ 934 try: 935 if kpt_line is None: 936 kpath = HighSymmKpath(structure) 937 kpt_line = [ 938 Kpoint(k, structure.lattice.reciprocal_lattice) 939 for k in kpath.get_kpoints(coords_are_cartesian=False)[0] 940 ] 941 labels_dict = {l: k for k, l in zip(*kpath.get_kpoints(coords_are_cartesian=False)) if l} 942 kpt_line = [kp.frac_coords for kp in kpt_line] 943 elif isinstance(kpt_line[0], Kpoint): 944 kpt_line = [kp.frac_coords for kp in kpt_line] 945 labels_dict = {k: labels_dict[k].frac_coords for k in labels_dict} 946 947 idx_list = [] 948 # kpt_dense=np.array([kp for kp in self._bz_kpoints]) 949 for i, kp in enumerate(kpt_line): 950 w = [] 951 prec = 1e-05 952 while len(w) == 0: 953 w = np.where(np.all(np.abs(kp - self._bz_kpoints) < [prec] * 3, axis=1))[0] 954 prec *= 10 955 956 # print( prec ) 957 idx_list.append([i, w[0]]) 958 959 # if len(w)>0: 960 # idx_list.append([i,w[0]]) 961 # else: 962 # w=np.where(np.all(np.abs(kp.frac_coords-self._bz_kpoints) 963 # <[1e-04,1e-04,1e-04],axis=1))[0] 964 # idx_list.append([i,w[0]]) 965 966 idx_list = np.array(idx_list) 967 # print( idx_list.shape ) 968 969 bands_dict = {Spin.up: (self._bz_bands * Energy(1, "Ry").to("eV") + efermi).T[:, idx_list[:, 1]].tolist()} 970 # bz_kpoints = bz_kpoints[idx_list[:,1]].tolist() 971 972 sbs = BandStructureSymmLine( 973 kpt_line, 974 bands_dict, 975 structure.lattice.reciprocal_lattice, 976 efermi, 977 labels_dict=labels_dict, 978 ) 979 980 return sbs 981 982 except Exception: 983 raise BoltztrapError( 984 "Bands are not in output of BoltzTraP.\nBolztrapRunner must " "be run with run_type=BANDS" 985 ) 986 987 @staticmethod 988 def check_acc_bzt_bands(sbs_bz, sbs_ref, warn_thr=(0.03, 0.03)): 989 """ 990 Compare sbs_bz BandStructureSymmLine calculated with boltztrap with 991 the sbs_ref BandStructureSymmLine as reference (from MP for 992 instance), computing correlation and energy difference for eight bands 993 around the gap (semiconductors) or fermi level (metals). 994 warn_thr is a threshold to get a warning in the accuracy of Boltztap 995 interpolated bands. 996 Return a dictionary with these keys: 997 - "N": the index of the band compared; inside each there are: 998 - "Corr": correlation coefficient for the 8 compared bands 999 - "Dist": energy distance for the 8 compared bands 1000 - "branch_name": energy distance for that branch 1001 - "avg_corr": average of correlation coefficient over the 8 bands 1002 - "avg_dist": average of energy distance over the 8 bands 1003 - "nb_list": list of indexes of the 8 compared bands 1004 - "acc_thr": list of two float corresponing to the two warning 1005 thresholds in input 1006 - "acc_err": list of two bools: 1007 True if the avg_corr > warn_thr[0], and 1008 True if the avg_dist > warn_thr[1] 1009 See also compare_sym_bands function doc 1010 """ 1011 if not sbs_ref.is_metal() and not sbs_bz.is_metal(): 1012 vbm_idx = sbs_bz.get_vbm()["band_index"][Spin.up][-1] 1013 cbm_idx = sbs_bz.get_cbm()["band_index"][Spin.up][0] 1014 nb_list = range(vbm_idx - 3, cbm_idx + 4) 1015 1016 else: 1017 bnd_around_efermi = [] 1018 delta = 0 1019 spin = list(sbs_bz.bands.keys())[0] 1020 while len(bnd_around_efermi) < 8 and delta < 100: 1021 delta += 0.1 1022 bnd_around_efermi = [] 1023 for nb in range(len(sbs_bz.bands[spin])): 1024 for kp in range(len(sbs_bz.bands[spin][nb])): 1025 if abs(sbs_bz.bands[spin][nb][kp] - sbs_bz.efermi) < delta: 1026 bnd_around_efermi.append(nb) 1027 break 1028 if len(bnd_around_efermi) < 8: 1029 print("Warning! check performed on " + str(len(bnd_around_efermi))) 1030 nb_list = bnd_around_efermi 1031 else: 1032 nb_list = bnd_around_efermi[:8] 1033 1034 # print(nb_list) 1035 bcheck = compare_sym_bands(sbs_bz, sbs_ref, nb_list) 1036 # print(bcheck) 1037 acc_err = [False, False] 1038 avg_corr = sum([item[1]["Corr"] for item in bcheck.items()]) / 8 1039 avg_distance = sum([item[1]["Dist"] for item in bcheck.items()]) / 8 1040 1041 if avg_corr > warn_thr[0]: 1042 acc_err[0] = True 1043 if avg_distance > warn_thr[0]: 1044 acc_err[1] = True 1045 1046 bcheck["avg_corr"] = avg_corr 1047 bcheck["avg_distance"] = avg_distance 1048 bcheck["acc_err"] = acc_err 1049 bcheck["acc_thr"] = warn_thr 1050 bcheck["nb_list"] = nb_list 1051 1052 if True in acc_err: 1053 print("Warning! some bands around gap are not accurate") 1054 1055 return bcheck 1056 1057 def get_seebeck(self, output="eigs", doping_levels=True): 1058 """ 1059 Gives the seebeck coefficient (microV/K) in either a 1060 full 3x3 tensor form, as 3 eigenvalues, or as the average value 1061 (trace/3.0) If doping_levels=True, the results are given at 1062 different p and n doping 1063 levels (given by self.doping), otherwise it is given as a series 1064 of electron chemical potential values 1065 1066 Args: 1067 output (string): the type of output. 'tensor' give the full 1068 3x3 tensor, 'eigs' its 3 eigenvalues and 1069 'average' the average of the three eigenvalues 1070 doping_levels (boolean): True for the results to be given at 1071 different doping levels, False for results 1072 at different electron chemical potentials 1073 1074 Returns: 1075 If doping_levels=True, a dictionary {temp:{'p':[],'n':[]}}. 1076 The 'p' links to Seebeck at p-type doping 1077 and 'n' to the Seebeck at n-type doping. Otherwise, returns a 1078 {temp:[]} dictionary 1079 The result contains either the sorted three eigenvalues of 1080 the symmetric 1081 Seebeck tensor (output='eigs') or a full tensor (3x3 array) ( 1082 output='tensor') or as an average 1083 (output='average'). 1084 1085 units are microV/K 1086 """ 1087 return BoltztrapAnalyzer._format_to_output(self._seebeck, self._seebeck_doping, output, doping_levels, 1e6) 1088 1089 def get_conductivity(self, output="eigs", doping_levels=True, relaxation_time=1e-14): 1090 """ 1091 Gives the conductivity (1/Ohm*m) in either a full 3x3 tensor 1092 form, as 3 eigenvalues, or as the average value 1093 (trace/3.0) If doping_levels=True, the results are given at 1094 different p and n doping 1095 levels (given by self.doping), otherwise it is given as a series 1096 of electron chemical potential values 1097 1098 Args: 1099 output (string): the type of output. 'tensor' give the full 1100 3x3 tensor, 'eigs' its 3 eigenvalues and 1101 'average' the average of the three eigenvalues 1102 doping_levels (boolean): True for the results to be given at 1103 different doping levels, False for results 1104 at different electron chemical potentials 1105 relaxation_time (float): constant relaxation time in secs 1106 1107 Returns: 1108 If doping_levels=True, a dictionary {temp:{'p':[],'n':[]}}. 1109 The 'p' links to conductivity 1110 at p-type doping and 'n' to the conductivity at n-type 1111 doping. Otherwise, 1112 returns a {temp:[]} dictionary. The result contains either 1113 the sorted three eigenvalues of the symmetric 1114 conductivity tensor (format='eigs') or a full tensor (3x3 1115 array) (output='tensor') or as an average 1116 (output='average'). 1117 The result includes a given constant relaxation time 1118 1119 units are 1/Ohm*m 1120 """ 1121 return BoltztrapAnalyzer._format_to_output( 1122 self._cond, self._cond_doping, output, doping_levels, relaxation_time 1123 ) 1124 1125 def get_power_factor(self, output="eigs", doping_levels=True, relaxation_time=1e-14): 1126 """ 1127 Gives the power factor (Seebeck^2 * conductivity) in units 1128 microW/(m*K^2) in either a full 3x3 tensor form, 1129 as 3 eigenvalues, or as the average value (trace/3.0) If 1130 doping_levels=True, the results are given at 1131 different p and n doping levels (given by self.doping), otherwise it 1132 is given as a series of 1133 electron chemical potential values 1134 1135 Args: 1136 output (string): the type of output. 'tensor' give the full 3x3 1137 tensor, 'eigs' its 3 eigenvalues and 1138 'average' the average of the three eigenvalues 1139 doping_levels (boolean): True for the results to be given at 1140 different doping levels, False for results 1141 at different electron chemical potentials 1142 relaxation_time (float): constant relaxation time in secs 1143 1144 Returns: 1145 If doping_levels=True, a dictionnary {temp:{'p':[],'n':[]}}. The 1146 'p' links to power factor 1147 at p-type doping and 'n' to the conductivity at n-type doping. 1148 Otherwise, 1149 returns a {temp:[]} dictionary. The result contains either the 1150 sorted three eigenvalues of the symmetric 1151 power factor tensor (format='eigs') or a full tensor (3x3 array) ( 1152 output='tensor') or as an average 1153 (output='average'). 1154 The result includes a given constant relaxation time 1155 1156 units are microW/(m K^2) 1157 """ 1158 result = None 1159 result_doping = None 1160 if doping_levels: 1161 result_doping = {doping: {t: [] for t in self._seebeck_doping[doping]} for doping in self._seebeck_doping} 1162 1163 for doping in result_doping: 1164 for t in result_doping[doping]: 1165 for i in range(len(self.doping[doping])): 1166 full_tensor = np.dot( 1167 self._cond_doping[doping][t][i], 1168 np.dot( 1169 self._seebeck_doping[doping][t][i], 1170 self._seebeck_doping[doping][t][i], 1171 ), 1172 ) 1173 result_doping[doping][t].append(full_tensor) 1174 1175 else: 1176 result = {t: [] for t in self._seebeck} 1177 for t in result: 1178 for i in range(len(self.mu_steps)): 1179 full_tensor = np.dot( 1180 self._cond[t][i], 1181 np.dot(self._seebeck[t][i], self._seebeck[t][i]), 1182 ) 1183 result[t].append(full_tensor) 1184 1185 return BoltztrapAnalyzer._format_to_output( 1186 result, result_doping, output, doping_levels, multi=1e6 * relaxation_time 1187 ) 1188 1189 def get_thermal_conductivity(self, output="eigs", doping_levels=True, k_el=True, relaxation_time=1e-14): 1190 """ 1191 Gives the electronic part of the thermal conductivity in either a 1192 full 3x3 tensor form, 1193 as 3 eigenvalues, or as the average value (trace/3.0) If 1194 doping_levels=True, the results are given at 1195 different p and n doping levels (given by self.doping), otherwise it 1196 is given as a series of 1197 electron chemical potential values 1198 1199 Args: 1200 output (string): the type of output. 'tensor' give the full 3x3 1201 tensor, 'eigs' its 3 eigenvalues and 1202 'average' the average of the three eigenvalues 1203 doping_levels (boolean): True for the results to be given at 1204 different doping levels, False for results 1205 at different electron chemical potentials 1206 k_el (boolean): True for k_0-PF*T, False for k_0 1207 relaxation_time (float): constant relaxation time in secs 1208 1209 Returns: 1210 If doping_levels=True, a dictionary {temp:{'p':[],'n':[]}}. The 1211 'p' links to thermal conductivity 1212 at p-type doping and 'n' to the thermal conductivity at n-type 1213 doping. Otherwise, 1214 returns a {temp:[]} dictionary. The result contains either the 1215 sorted three eigenvalues of the symmetric 1216 conductivity tensor (format='eigs') or a full tensor (3x3 array) ( 1217 output='tensor') or as an average 1218 (output='average'). 1219 The result includes a given constant relaxation time 1220 1221 units are W/mK 1222 """ 1223 result = None 1224 result_doping = None 1225 if doping_levels: 1226 result_doping = {doping: {t: [] for t in self._seebeck_doping[doping]} for doping in self._seebeck_doping} 1227 for doping in result_doping: 1228 for t in result_doping[doping]: 1229 for i in range(len(self.doping[doping])): 1230 if k_el: 1231 pf_tensor = np.dot( 1232 self._cond_doping[doping][t][i], 1233 np.dot( 1234 self._seebeck_doping[doping][t][i], 1235 self._seebeck_doping[doping][t][i], 1236 ), 1237 ) 1238 result_doping[doping][t].append((self._kappa_doping[doping][t][i] - pf_tensor * t)) 1239 else: 1240 result_doping[doping][t].append((self._kappa_doping[doping][t][i])) 1241 else: 1242 result = {t: [] for t in self._seebeck} 1243 for t in result: 1244 for i in range(len(self.mu_steps)): 1245 if k_el: 1246 pf_tensor = np.dot( 1247 self._cond[t][i], 1248 np.dot(self._seebeck[t][i], self._seebeck[t][i]), 1249 ) 1250 result[t].append((self._kappa[t][i] - pf_tensor * t)) 1251 else: 1252 result[t].append((self._kappa[t][i])) 1253 1254 return BoltztrapAnalyzer._format_to_output(result, result_doping, output, doping_levels, multi=relaxation_time) 1255 1256 def get_zt(self, output="eigs", doping_levels=True, relaxation_time=1e-14, kl=1.0): 1257 """ 1258 Gives the ZT coefficient (S^2*cond*T/thermal cond) in either a full 1259 3x3 tensor form, 1260 as 3 eigenvalues, or as the average value (trace/3.0) If 1261 doping_levels=True, the results are given at 1262 different p and n doping levels (given by self.doping), otherwise it 1263 is given as a series of 1264 electron chemical potential values. We assume a constant relaxation 1265 time and a constant 1266 lattice thermal conductivity 1267 1268 Args: 1269 output (string): the type of output. 'tensor' give the full 3x3 1270 tensor, 'eigs' its 3 eigenvalues and 1271 'average' the average of the three eigenvalues 1272 doping_levels (boolean): True for the results to be given at 1273 different doping levels, False for results 1274 at different electron chemical potentials 1275 relaxation_time (float): constant relaxation time in secs 1276 k_l (float): lattice thermal cond in W/(m*K) 1277 1278 Returns: 1279 If doping_levels=True, a dictionary {temp:{'p':[],'n':[]}}. The 1280 'p' links to ZT 1281 at p-type doping and 'n' to the ZT at n-type doping. Otherwise, 1282 returns a {temp:[]} dictionary. The result contains either the 1283 sorted three eigenvalues of the symmetric 1284 ZT tensor (format='eigs') or a full tensor (3x3 array) ( 1285 output='tensor') or as an average 1286 (output='average'). 1287 The result includes a given constant relaxation time and lattice 1288 thermal conductivity 1289 """ 1290 result = None 1291 result_doping = None 1292 if doping_levels: 1293 result_doping = {doping: {t: [] for t in self._seebeck_doping[doping]} for doping in self._seebeck_doping} 1294 1295 for doping in result_doping: 1296 for t in result_doping[doping]: 1297 for i in range(len(self.doping[doping])): 1298 pf_tensor = np.dot( 1299 self._cond_doping[doping][t][i], 1300 np.dot( 1301 self._seebeck_doping[doping][t][i], 1302 self._seebeck_doping[doping][t][i], 1303 ), 1304 ) 1305 thermal_conduct = (self._kappa_doping[doping][t][i] - pf_tensor * t) * relaxation_time 1306 result_doping[doping][t].append( 1307 np.dot( 1308 pf_tensor * relaxation_time * t, 1309 np.linalg.inv(thermal_conduct + kl * np.eye(3, 3)), 1310 ) 1311 ) 1312 else: 1313 result = {t: [] for t in self._seebeck} 1314 for t in result: 1315 for i in range(len(self.mu_steps)): 1316 pf_tensor = np.dot( 1317 self._cond[t][i], 1318 np.dot(self._seebeck[t][i], self._seebeck[t][i]), 1319 ) 1320 thermal_conduct = (self._kappa[t][i] - pf_tensor * t) * relaxation_time 1321 result[t].append( 1322 np.dot( 1323 pf_tensor * relaxation_time * t, 1324 np.linalg.inv(thermal_conduct + kl * np.eye(3, 3)), 1325 ) 1326 ) 1327 1328 return BoltztrapAnalyzer._format_to_output(result, result_doping, output, doping_levels) 1329 1330 def get_average_eff_mass(self, output="eigs", doping_levels=True): 1331 """ 1332 Gives the average effective mass tensor. We call it average because 1333 it takes into account all the bands 1334 and regions in the Brillouin zone. This is different than the standard 1335 textbook effective mass which relates 1336 often to only one (parabolic) band. 1337 The average effective mass tensor is defined as the integrated 1338 average of the second derivative of E(k) 1339 This effective mass tensor takes into account: 1340 -non-parabolicity 1341 -multiple extrema 1342 -multiple bands 1343 1344 For more information about it. See: 1345 1346 Hautier, G., Miglio, A., Waroquiers, D., Rignanese, G., & Gonze, 1347 X. (2014). 1348 How Does Chemistry Influence Electron Effective Mass in Oxides? 1349 A High-Throughput Computational Analysis. Chemistry of Materials, 1350 26(19), 5447–5458. doi:10.1021/cm404079a 1351 1352 or 1353 1354 Hautier, G., Miglio, A., Ceder, G., Rignanese, G.-M., & Gonze, 1355 X. (2013). 1356 Identification and design principles of low hole effective mass 1357 p-type transparent conducting oxides. 1358 Nature Communications, 4, 2292. doi:10.1038/ncomms3292 1359 1360 Depending on the value of output, we have either the full 3x3 1361 effective mass tensor, 1362 its 3 eigenvalues or an average 1363 1364 Args: 1365 output (string): 'eigs' for eigenvalues, 'tensor' for the full 1366 tensor and 'average' for an average (trace/3) 1367 doping_levels (boolean): True for the results to be given at 1368 different doping levels, False for results 1369 at different electron chemical potentials 1370 Returns: 1371 If doping_levels=True,a dictionary {'p':{temp:[]},'n':{temp:[]}} 1372 with an array of effective mass tensor, eigenvalues of average 1373 value (depending on output) for each temperature and for each 1374 doping level. 1375 The 'p' links to hole effective mass tensor and 'n' to electron 1376 effective mass tensor. 1377 """ 1378 result = None 1379 result_doping = None 1380 conc = self.get_carrier_concentration() 1381 if doping_levels: 1382 result_doping = {doping: {t: [] for t in self._cond_doping[doping]} for doping in self.doping} 1383 for doping in result_doping: 1384 for temp in result_doping[doping]: 1385 for i in range(len(self.doping[doping])): 1386 try: 1387 result_doping[doping][temp].append( 1388 np.linalg.inv(np.array(self._cond_doping[doping][temp][i])) 1389 * self.doping[doping][i] 1390 * 10 ** 6 1391 * constants.e ** 2 1392 / constants.m_e 1393 ) 1394 except np.linalg.LinAlgError: 1395 pass 1396 else: 1397 result = {t: [] for t in self._seebeck} 1398 for temp in result: 1399 for i in range(len(self.mu_steps)): 1400 try: 1401 cond_inv = np.linalg.inv(np.array(self._cond[temp][i])) 1402 except np.linalg.LinAlgError: 1403 pass 1404 result[temp].append(cond_inv * conc[temp][i] * 10 ** 6 * constants.e ** 2 / constants.m_e) 1405 1406 return BoltztrapAnalyzer._format_to_output(result, result_doping, output, doping_levels) 1407 1408 def get_seebeck_eff_mass(self, output="average", temp=300, doping_levels=False, Lambda=0.5): 1409 """ 1410 Seebeck effective mass calculated as explained in Ref. 1411 Gibbs, Z. M. et al., Effective mass and fermi surface complexity factor 1412 from ab initio band structure calculations. 1413 npj Computational Materials 3, 8 (2017). 1414 1415 Args: 1416 output: 'average' returns the seebeck effective mass calculated using 1417 the average of the three diagonal components of the seebeck tensor. 1418 'tensor' returns the seebeck effective mass respect to the three 1419 diagonal components of the seebeck tensor. 1420 doping_levels: False means that the seebeck effective mass is calculated 1421 for every value of the chemical potential 1422 True means that the seebeck effective mass is calculated 1423 for every value of the doping levels for both n and p types 1424 temp: temperature of calculated seebeck. 1425 Lambda: fitting parameter used to model the scattering (0.5 means constant 1426 relaxation time). 1427 Returns: 1428 a list of values for the seebeck effective mass w.r.t the chemical potential, 1429 if doping_levels is set at False; 1430 a dict with n an p keys that contain a list of values for the seebeck effective 1431 mass w.r.t the doping levels, if doping_levels is set at True; 1432 if 'tensor' is selected, each element of the lists is a list containing 1433 the three components of the seebeck effective mass. 1434 """ 1435 1436 if doping_levels: 1437 sbk_mass = {} 1438 for dt in ("n", "p"): 1439 conc = self.doping[dt] 1440 seebeck = self.get_seebeck(output=output, doping_levels=True)[dt][temp] 1441 sbk_mass[dt] = [] 1442 for i, c in enumerate(conc): 1443 if output == "average": 1444 sbk_mass[dt].append(seebeck_eff_mass_from_seebeck_carr(abs(seebeck[i]), c, temp, Lambda)) 1445 elif output == "tensor": 1446 sbk_mass[dt].append([]) 1447 for j in range(3): 1448 sbk_mass[dt][-1].append( 1449 seebeck_eff_mass_from_seebeck_carr(abs(seebeck[i][j][j]), c, temp, Lambda) 1450 ) 1451 1452 else: 1453 seebeck = self.get_seebeck(output=output, doping_levels=False)[temp] 1454 conc = self.get_carrier_concentration()[temp] 1455 sbk_mass = [] 1456 for i, c in enumerate(conc): 1457 if output == "average": 1458 sbk_mass.append(seebeck_eff_mass_from_seebeck_carr(abs(seebeck[i]), c, temp, Lambda)) 1459 elif output == "tensor": 1460 sbk_mass.append([]) 1461 for j in range(3): 1462 sbk_mass[-1].append(seebeck_eff_mass_from_seebeck_carr(abs(seebeck[i][j][j]), c, temp, Lambda)) 1463 return sbk_mass 1464 1465 def get_complexity_factor(self, output="average", temp=300, doping_levels=False, Lambda=0.5): 1466 """ 1467 Fermi surface complexity factor respect to calculated as explained in Ref. 1468 Gibbs, Z. M. et al., Effective mass and fermi surface complexity factor 1469 from ab initio band structure calculations. 1470 npj Computational Materials 3, 8 (2017). 1471 1472 Args: 1473 output: 'average' returns the complexity factor calculated using the average 1474 of the three diagonal components of the seebeck and conductivity tensors. 1475 'tensor' returns the complexity factor respect to the three 1476 diagonal components of seebeck and conductivity tensors. 1477 doping_levels: False means that the complexity factor is calculated 1478 for every value of the chemical potential 1479 True means that the complexity factor is calculated 1480 for every value of the doping levels for both n and p types 1481 temp: temperature of calculated seebeck and conductivity. 1482 Lambda: fitting parameter used to model the scattering (0.5 means constant 1483 relaxation time). 1484 Returns: 1485 a list of values for the complexity factor w.r.t the chemical potential, 1486 if doping_levels is set at False; 1487 a dict with n an p keys that contain a list of values for the complexity factor 1488 w.r.t the doping levels, if doping_levels is set at True; 1489 if 'tensor' is selected, each element of the lists is a list containing 1490 the three components of the complexity factor. 1491 """ 1492 1493 if doping_levels: 1494 cmplx_fact = {} 1495 for dt in ("n", "p"): 1496 sbk_mass = self.get_seebeck_eff_mass(output, temp, True, Lambda)[dt] 1497 cond_mass = self.get_average_eff_mass(output=output, doping_levels=True)[dt][temp] 1498 1499 if output == "average": 1500 cmplx_fact[dt] = [(m_s / abs(m_c)) ** 1.5 for m_s, m_c in zip(sbk_mass, cond_mass)] 1501 elif output == "tensor": 1502 cmplx_fact[dt] = [] 1503 for i, sm in enumerate(sbk_mass): 1504 cmplx_fact[dt].append([]) 1505 for j in range(3): 1506 cmplx_fact[dt][-1].append((sm[j] / abs(cond_mass[i][j][j])) ** 1.5) 1507 1508 else: 1509 sbk_mass = self.get_seebeck_eff_mass(output, temp, False, Lambda) 1510 cond_mass = self.get_average_eff_mass(output=output, doping_levels=False)[temp] 1511 1512 if output == "average": 1513 cmplx_fact = [(m_s / abs(m_c)) ** 1.5 for m_s, m_c in zip(sbk_mass, cond_mass)] 1514 elif output == "tensor": 1515 cmplx_fact = [] 1516 for i, sm in enumerate(sbk_mass): 1517 cmplx_fact.append([]) 1518 for j in range(3): 1519 cmplx_fact[-1].append((sm[j] / abs(cond_mass[i][j][j])) ** 1.5) 1520 1521 return cmplx_fact 1522 1523 def get_extreme( 1524 self, 1525 target_prop, 1526 maximize=True, 1527 min_temp=None, 1528 max_temp=None, 1529 min_doping=None, 1530 max_doping=None, 1531 isotropy_tolerance=0.05, 1532 use_average=True, 1533 ): 1534 """ 1535 This method takes in eigenvalues over a range of carriers, 1536 temperatures, and doping levels, and tells you what is the "best" 1537 value that can be achieved for the given target_property. Note that 1538 this method searches the doping dict only, not the full mu dict. 1539 1540 Args: 1541 target_prop: target property, i.e. "seebeck", "power factor", 1542 "conductivity", "kappa", or "zt" 1543 maximize: True to maximize, False to minimize (e.g. kappa) 1544 min_temp: minimum temperature allowed 1545 max_temp: maximum temperature allowed 1546 min_doping: minimum doping allowed (e.g., 1E18) 1547 max_doping: maximum doping allowed (e.g., 1E20) 1548 isotropy_tolerance: tolerance for isotropic (0.05 = 5%) 1549 use_average: True for avg of eigenval, False for max eigenval 1550 1551 Returns: 1552 A dictionary with keys {"p", "n", "best"} with sub-keys: 1553 {"value", "temperature", "doping", "isotropic"} 1554 1555 """ 1556 1557 def is_isotropic(x, isotropy_tolerance): 1558 """ 1559 Internal method to tell you if 3-vector "x" is isotropic 1560 1561 Args: 1562 x: the vector to determine isotropy for 1563 isotropy_tolerance: tolerance, e.g. 0.05 is 5% 1564 """ 1565 if len(x) != 3: 1566 raise ValueError("Invalid input to is_isotropic!") 1567 1568 st = sorted(x) 1569 return bool( 1570 all([st[0], st[1], st[2]]) 1571 and (abs((st[1] - st[0]) / st[1]) <= isotropy_tolerance) 1572 and (abs((st[2] - st[0])) / st[2] <= isotropy_tolerance) 1573 and (abs((st[2] - st[1]) / st[2]) <= isotropy_tolerance) 1574 ) 1575 1576 if target_prop.lower() == "seebeck": 1577 d = self.get_seebeck(output="eigs", doping_levels=True) 1578 1579 elif target_prop.lower() == "power factor": 1580 d = self.get_power_factor(output="eigs", doping_levels=True) 1581 1582 elif target_prop.lower() == "conductivity": 1583 d = self.get_conductivity(output="eigs", doping_levels=True) 1584 1585 elif target_prop.lower() == "kappa": 1586 d = self.get_thermal_conductivity(output="eigs", doping_levels=True) 1587 elif target_prop.lower() == "zt": 1588 d = self.get_zt(output="eigs", doping_levels=True) 1589 1590 else: 1591 raise ValueError("Target property: {} not recognized!".format(target_prop)) 1592 1593 absval = True # take the absolute value of properties 1594 1595 x_val = None 1596 x_temp = None 1597 x_doping = None 1598 x_isotropic = None 1599 output = {} 1600 1601 min_temp = min_temp or 0 1602 max_temp = max_temp or float("inf") 1603 min_doping = min_doping or 0 1604 max_doping = max_doping or float("inf") 1605 1606 for pn in ("p", "n"): 1607 for t in d[pn]: # temperatures 1608 if min_temp <= float(t) <= max_temp: 1609 for didx, evs in enumerate(d[pn][t]): 1610 doping_lvl = self.doping[pn][didx] 1611 if min_doping <= doping_lvl <= max_doping: 1612 isotropic = is_isotropic(evs, isotropy_tolerance) 1613 if absval: 1614 evs = [abs(x) for x in evs] 1615 if use_average: 1616 val = float(sum(evs)) / len(evs) 1617 else: 1618 val = max(evs) 1619 if x_val is None or (val > x_val and maximize) or (val < x_val and not maximize): 1620 x_val = val 1621 x_temp = t 1622 x_doping = doping_lvl 1623 x_isotropic = isotropic 1624 1625 output[pn] = { 1626 "value": x_val, 1627 "temperature": x_temp, 1628 "doping": x_doping, 1629 "isotropic": x_isotropic, 1630 } 1631 x_val = None 1632 1633 if maximize: 1634 max_type = "p" if output["p"]["value"] >= output["n"]["value"] else "n" 1635 else: 1636 max_type = "p" if output["p"]["value"] <= output["n"]["value"] else "n" 1637 1638 output["best"] = output[max_type] 1639 output["best"]["carrier_type"] = max_type 1640 1641 return output 1642 1643 @staticmethod 1644 def _format_to_output(tensor, tensor_doping, output, doping_levels, multi=1.0): 1645 if doping_levels: 1646 full_tensor = tensor_doping 1647 result = {doping: {t: [] for t in tensor_doping[doping]} for doping in tensor_doping} 1648 for doping in full_tensor: 1649 for temp in full_tensor[doping]: 1650 for i in range(len(full_tensor[doping][temp])): 1651 if output in ["eig", "eigs"]: 1652 result[doping][temp].append(sorted(np.linalg.eigh(full_tensor[doping][temp][i])[0] * multi)) 1653 elif output == "tensor": 1654 result[doping][temp].append(np.array(full_tensor[doping][temp][i]) * multi) 1655 elif output == "average": 1656 result[doping][temp].append( 1657 ( 1658 full_tensor[doping][temp][i][0][0] 1659 + full_tensor[doping][temp][i][1][1] 1660 + full_tensor[doping][temp][i][2][2] 1661 ) 1662 * multi 1663 / 3.0 1664 ) 1665 else: 1666 raise ValueError("Unknown output format: " "{}".format(output)) 1667 else: 1668 full_tensor = tensor 1669 result = {t: [] for t in tensor} 1670 for temp in full_tensor: 1671 for i in range(len(tensor[temp])): 1672 if output in ["eig", "eigs"]: 1673 result[temp].append(sorted(np.linalg.eigh(full_tensor[temp][i])[0] * multi)) 1674 elif output == "tensor": 1675 result[temp].append(np.array(full_tensor[temp][i]) * multi) 1676 elif output == "average": 1677 result[temp].append( 1678 (full_tensor[temp][i][0][0] + full_tensor[temp][i][1][1] + full_tensor[temp][i][2][2]) 1679 * multi 1680 / 3.0 1681 ) 1682 else: 1683 raise ValueError("Unknown output format: {}".format(output)) 1684 return result 1685 1686 def get_complete_dos(self, structure, analyzer_for_second_spin=None): 1687 """ 1688 Gives a CompleteDos object with the DOS from the interpolated 1689 projected band structure 1690 Args: 1691 the structure (necessary to identify sites for projection) 1692 analyzer_for_second_spin must be specified to have a 1693 CompleteDos with both Spin components 1694 Returns: 1695 a CompleteDos object 1696 Example of use in case of spin polarized case: 1697 1698 BoltztrapRunner(bs=bs,nelec=10,run_type="DOS",spin=1).run(path_dir='dos_up/') 1699 an_up=BoltztrapAnalyzer.from_files("dos_up/boltztrap/",dos_spin=1) 1700 1701 BoltztrapRunner(bs=bs,nelec=10,run_type="DOS",spin=-1).run(path_dir='dos_dw/') 1702 an_dw=BoltztrapAnalyzer.from_files("dos_dw/boltztrap/",dos_spin=-1) 1703 1704 cdos=an_up.get_complete_dos(bs.structure,an_dw) 1705 1706 """ 1707 pdoss = {} 1708 spin_1 = list(self.dos.densities.keys())[0] 1709 1710 if analyzer_for_second_spin: 1711 if not np.all(self.dos.energies == analyzer_for_second_spin.dos.energies): 1712 raise BoltztrapError("Dos merging error: energies of the two dos are different") 1713 1714 spin_2 = list(analyzer_for_second_spin.dos.densities.keys())[0] 1715 if spin_1 == spin_2: 1716 raise BoltztrapError("Dos merging error: spin component are the same") 1717 1718 for s in self._dos_partial: 1719 if structure.sites[int(s)] not in pdoss: 1720 pdoss[structure.sites[int(s)]] = {} 1721 for o in self._dos_partial[s]: 1722 if Orbital[o] not in pdoss[structure.sites[int(s)]]: 1723 pdoss[structure.sites[int(s)]][Orbital[o]] = {} 1724 pdoss[structure.sites[int(s)]][Orbital[o]][spin_1] = self._dos_partial[s][o] 1725 if analyzer_for_second_spin: 1726 pdoss[structure.sites[int(s)]][Orbital[o]][spin_2] = analyzer_for_second_spin._dos_partial[s][o] 1727 if analyzer_for_second_spin: 1728 tdos = Dos( 1729 self.dos.efermi, 1730 self.dos.energies, 1731 { 1732 spin_1: self.dos.densities[spin_1], 1733 spin_2: analyzer_for_second_spin.dos.densities[spin_2], 1734 }, 1735 ) 1736 else: 1737 tdos = self.dos 1738 1739 return CompleteDos(structure, total_dos=tdos, pdoss=pdoss) 1740 1741 def get_mu_bounds(self, temp=300): 1742 """ 1743 :param temp: Temperature. 1744 :return: The chemical potential bounds at that temperature. 1745 """ 1746 return min(self.mu_doping["p"][temp]), max(self.mu_doping["n"][temp]) 1747 1748 def get_carrier_concentration(self): 1749 """ 1750 gives the carrier concentration (in cm^-3) 1751 1752 Returns 1753 a dictionary {temp:[]} with an array of carrier concentration 1754 (in cm^-3) at each temperature 1755 The array relates to each step of electron chemical potential 1756 """ 1757 1758 return {temp: [1e24 * i / self.vol for i in self._carrier_conc[temp]] for temp in self._carrier_conc} 1759 1760 def get_hall_carrier_concentration(self): 1761 """ 1762 gives the Hall carrier concentration (in cm^-3). This is the trace of 1763 the Hall tensor (see Boltztrap source code) Hall carrier concentration 1764 are not always exactly the same than carrier concentration. 1765 1766 Returns 1767 a dictionary {temp:[]} with an array of Hall carrier concentration 1768 (in cm^-3) at each temperature The array relates to each step of 1769 electron chemical potential 1770 """ 1771 result = {temp: [] for temp in self._hall} 1772 for temp in self._hall: 1773 for i in self._hall[temp]: 1774 trace = (i[1][2][0] + i[2][0][1] + i[0][1][2]) / 3.0 1775 if trace != 0.0: 1776 result[temp].append(1e-6 / (trace * constants.e)) 1777 else: 1778 result[temp].append(0.0) 1779 return result 1780 1781 @staticmethod 1782 def parse_outputtrans(path_dir): 1783 """ 1784 Parses .outputtrans file 1785 1786 Args: 1787 path_dir: dir containing boltztrap.outputtrans 1788 1789 Returns: 1790 tuple - (run_type, warning, efermi, gap, doping_levels) 1791 1792 """ 1793 run_type = None 1794 warning = None 1795 efermi = None 1796 gap = None 1797 doping_levels = [] 1798 1799 with open(os.path.join(path_dir, "boltztrap.outputtrans"), "r") as f: 1800 for line in f: 1801 if "WARNING" in line: 1802 warning = line 1803 elif "Calc type:" in line: 1804 run_type = line.split()[-1] 1805 elif line.startswith("VBM"): 1806 efermi = Energy(line.split()[1], "Ry").to("eV") 1807 elif line.startswith("Egap:"): 1808 gap = Energy(float(line.split()[1]), "Ry").to("eV") 1809 elif line.startswith("Doping level number"): 1810 doping_levels.append(float(line.split()[6])) 1811 1812 return run_type, warning, efermi, gap, doping_levels 1813 1814 @staticmethod 1815 def parse_transdos(path_dir, efermi, dos_spin=1, trim_dos=False): 1816 """ 1817 Parses .transdos (total DOS) and .transdos_x_y (partial DOS) files 1818 Args: 1819 path_dir: (str) dir containing DOS files 1820 efermi: (float) Fermi energy 1821 dos_spin: (int) -1 for spin down, +1 for spin up 1822 trim_dos: (bool) whether to post-process / trim DOS 1823 1824 Returns: 1825 tuple - (DOS, dict of partial DOS) 1826 """ 1827 1828 data_dos = {"total": [], "partial": {}} 1829 # parse the total DOS data 1830 # format is energy, DOS, integrated DOS 1831 with open(os.path.join(path_dir, "boltztrap.transdos"), "r") as f: 1832 count_series = 0 # TODO: why is count_series needed? 1833 for line in f: 1834 if line.lstrip().startswith("#"): 1835 count_series += 1 1836 if count_series > 1: 1837 break 1838 else: 1839 data_dos["total"].append( 1840 [ 1841 Energy(float(line.split()[0]), "Ry").to("eV"), 1842 float(line.split()[1]), 1843 ] 1844 ) 1845 1846 lw_l = 0 1847 hg_l = -len(data_dos["total"]) 1848 if trim_dos: 1849 # Francesco knows what this does 1850 # It has something to do with a trick of adding fake energies 1851 # at the endpoints of the DOS, and then re-trimming it. This is 1852 # to get the same energy scale for up and down spin DOS. 1853 tmp_data = np.array(data_dos["total"]) 1854 tmp_den = np.trim_zeros(tmp_data[:, 1], "f")[1:] 1855 lw_l = len(tmp_data[:, 1]) - len(tmp_den) 1856 tmp_ene = tmp_data[lw_l:, 0] 1857 tmp_den = np.trim_zeros(tmp_den, "b")[:-1] 1858 hg_l = len(tmp_ene) - len(tmp_den) 1859 tmp_ene = tmp_ene[:-hg_l] 1860 tmp_data = np.vstack((tmp_ene, tmp_den)).T 1861 data_dos["total"] = tmp_data.tolist() 1862 1863 # parse partial DOS data 1864 for file_name in os.listdir(path_dir): 1865 if file_name.endswith("transdos") and file_name != "boltztrap.transdos": 1866 tokens = file_name.split(".")[1].split("_") 1867 site = tokens[1] 1868 orb = "_".join(tokens[2:]) 1869 with open(os.path.join(path_dir, file_name), "r") as f: 1870 for line in f: 1871 if not line.lstrip().startswith(" #"): 1872 if site not in data_dos["partial"]: 1873 data_dos["partial"][site] = {} 1874 if orb not in data_dos["partial"][site]: 1875 data_dos["partial"][site][orb] = [] 1876 data_dos["partial"][site][orb].append(float(line.split()[1])) 1877 data_dos["partial"][site][orb] = data_dos["partial"][site][orb][lw_l:-hg_l] 1878 1879 dos_full = {"energy": [], "density": []} 1880 1881 for t in data_dos["total"]: 1882 dos_full["energy"].append(t[0]) 1883 dos_full["density"].append(t[1]) 1884 1885 dos = Dos(efermi, dos_full["energy"], {Spin(dos_spin): dos_full["density"]}) 1886 dos_partial = data_dos["partial"] # TODO: make this real DOS object? 1887 1888 return dos, dos_partial 1889 1890 @staticmethod 1891 def parse_intrans(path_dir): 1892 """ 1893 Parses boltztrap.intrans mainly to extract the value of scissor applied 1894 to the bands or some other inputs 1895 1896 Args: 1897 path_dir: (str) dir containing the boltztrap.intrans file 1898 1899 Returns: 1900 intrans (dict): a dictionary containing various inputs that had 1901 been used in the Boltztrap run. 1902 """ 1903 intrans = {} 1904 with open(os.path.join(path_dir, "boltztrap.intrans"), "r") as f: 1905 for line in f: 1906 if "iskip" in line: 1907 intrans["scissor"] = Energy(float(line.split(" ")[3]), "Ry").to("eV") 1908 if "HISTO" in line or "TETRA" in line: 1909 intrans["dos_type"] = line[:-1] 1910 return intrans 1911 1912 @staticmethod 1913 def parse_struct(path_dir): 1914 """ 1915 Parses boltztrap.struct file (only the volume) 1916 Args: 1917 path_dir: (str) dir containing the boltztrap.struct file 1918 1919 Returns: 1920 (float) volume 1921 """ 1922 with open(os.path.join(path_dir, "boltztrap.struct"), "r") as f: 1923 tokens = f.readlines() 1924 return Lattice( 1925 [[Length(float(tokens[i].split()[j]), "bohr").to("ang") for j in range(3)] for i in range(1, 4)] 1926 ).volume 1927 1928 @staticmethod 1929 def parse_cond_and_hall(path_dir, doping_levels=None): 1930 """ 1931 Parses the conductivity and Hall tensors 1932 Args: 1933 path_dir: Path containing .condtens / .halltens files 1934 doping_levels: ([float]) - doping lvls, parse outtrans to get this 1935 1936 Returns: 1937 mu_steps, cond, seebeck, kappa, hall, pn_doping_levels, 1938 mu_doping, seebeck_doping, cond_doping, kappa_doping, 1939 hall_doping, carrier_conc 1940 """ 1941 1942 # Step 1: parse raw data but do not convert to final format 1943 t_steps = set() 1944 mu_steps = set() 1945 data_full = [] 1946 data_hall = [] 1947 data_doping_full = [] 1948 data_doping_hall = [] 1949 doping_levels = doping_levels or [] 1950 1951 # parse the full conductivity/Seebeck/kappa0/etc data 1952 # also initialize t_steps and mu_steps 1953 with open(os.path.join(path_dir, "boltztrap.condtens"), "r") as f: 1954 for line in f: 1955 if not line.startswith("#"): 1956 mu_steps.add(float(line.split()[0])) 1957 t_steps.add(int(float(line.split()[1]))) 1958 data_full.append([float(c) for c in line.split()]) 1959 1960 # parse the full Hall tensor 1961 with open(os.path.join(path_dir, "boltztrap.halltens"), "r") as f: 1962 for line in f: 1963 if not line.startswith("#"): 1964 data_hall.append([float(c) for c in line.split()]) 1965 1966 if len(doping_levels) != 0: 1967 # parse doping levels version of full cond. tensor, etc. 1968 with open(os.path.join(path_dir, "boltztrap.condtens_fixdoping"), "r") as f: 1969 for line in f: 1970 if not line.startswith("#") and len(line) > 2: 1971 data_doping_full.append([float(c) for c in line.split()]) 1972 1973 # parse doping levels version of full hall tensor 1974 with open(os.path.join(path_dir, "boltztrap.halltens_fixdoping"), "r") as f: 1975 for line in f: 1976 if not line.startswith("#") and len(line) > 2: 1977 data_doping_hall.append([float(c) for c in line.split()]) 1978 1979 # Step 2: convert raw data to final format 1980 1981 # sort t and mu_steps (b/c they are sets not lists) 1982 # and convert to correct energy 1983 t_steps = sorted(t_steps) 1984 mu_steps = sorted([Energy(m, "Ry").to("eV") for m in mu_steps]) 1985 1986 # initialize output variables - could use defaultdict instead 1987 # I am leaving things like this for clarity 1988 cond = {t: [] for t in t_steps} 1989 seebeck = {t: [] for t in t_steps} 1990 kappa = {t: [] for t in t_steps} 1991 hall = {t: [] for t in t_steps} 1992 carrier_conc = {t: [] for t in t_steps} 1993 1994 mu_doping = {"p": {t: [] for t in t_steps}, "n": {t: [] for t in t_steps}} 1995 seebeck_doping = {"p": {t: [] for t in t_steps}, "n": {t: [] for t in t_steps}} 1996 cond_doping = {"p": {t: [] for t in t_steps}, "n": {t: [] for t in t_steps}} 1997 kappa_doping = {"p": {t: [] for t in t_steps}, "n": {t: [] for t in t_steps}} 1998 hall_doping = {"p": {t: [] for t in t_steps}, "n": {t: [] for t in t_steps}} 1999 2000 # process doping levels 2001 pn_doping_levels = {"p": [], "n": []} 2002 for d in doping_levels: 2003 if d > 0: 2004 pn_doping_levels["p"].append(d) 2005 else: 2006 pn_doping_levels["n"].append(-d) 2007 2008 # process raw conductivity data, etc. 2009 for d in data_full: 2010 temp, doping = d[1], d[2] 2011 carrier_conc[temp].append(doping) 2012 2013 cond[temp].append(np.reshape(d[3:12], (3, 3)).tolist()) 2014 seebeck[temp].append(np.reshape(d[12:21], (3, 3)).tolist()) 2015 kappa[temp].append(np.reshape(d[21:30], (3, 3)).tolist()) 2016 2017 # process raw Hall data 2018 for d in data_hall: 2019 temp, doping = d[1], d[2] 2020 hall_tens = [ 2021 np.reshape(d[3:12], (3, 3)).tolist(), 2022 np.reshape(d[12:21], (3, 3)).tolist(), 2023 np.reshape(d[21:30], (3, 3)).tolist(), 2024 ] 2025 hall[temp].append(hall_tens) 2026 2027 # process doping conductivity data, etc. 2028 for d in data_doping_full: 2029 temp, doping, mu = d[0], d[1], d[-1] 2030 pn = "p" if doping > 0 else "n" 2031 mu_doping[pn][temp].append(Energy(mu, "Ry").to("eV")) 2032 cond_doping[pn][temp].append(np.reshape(d[2:11], (3, 3)).tolist()) 2033 seebeck_doping[pn][temp].append(np.reshape(d[11:20], (3, 3)).tolist()) 2034 kappa_doping[pn][temp].append(np.reshape(d[20:29], (3, 3)).tolist()) 2035 2036 # process doping Hall data 2037 for d in data_doping_hall: 2038 temp, doping, mu = d[0], d[1], d[-1] 2039 pn = "p" if doping > 0 else "n" 2040 hall_tens = [ 2041 np.reshape(d[2:11], (3, 3)).tolist(), 2042 np.reshape(d[11:20], (3, 3)).tolist(), 2043 np.reshape(d[20:29], (3, 3)).tolist(), 2044 ] 2045 hall_doping[pn][temp].append(hall_tens) 2046 2047 return ( 2048 mu_steps, 2049 cond, 2050 seebeck, 2051 kappa, 2052 hall, 2053 pn_doping_levels, 2054 mu_doping, 2055 seebeck_doping, 2056 cond_doping, 2057 kappa_doping, 2058 hall_doping, 2059 carrier_conc, 2060 ) 2061 2062 @staticmethod 2063 def from_files(path_dir, dos_spin=1): 2064 """ 2065 get a BoltztrapAnalyzer object from a set of files 2066 2067 Args: 2068 path_dir: directory where the boltztrap files are 2069 dos_spin: in DOS mode, set to 1 for spin up and -1 for spin down 2070 2071 Returns: 2072 a BoltztrapAnalyzer object 2073 2074 """ 2075 ( 2076 run_type, 2077 warning, 2078 efermi, 2079 gap, 2080 doping_levels, 2081 ) = BoltztrapAnalyzer.parse_outputtrans(path_dir) 2082 2083 vol = BoltztrapAnalyzer.parse_struct(path_dir) 2084 2085 intrans = BoltztrapAnalyzer.parse_intrans(path_dir) 2086 2087 if run_type == "BOLTZ": 2088 dos, pdos = BoltztrapAnalyzer.parse_transdos(path_dir, efermi, dos_spin=dos_spin, trim_dos=False) 2089 2090 ( 2091 mu_steps, 2092 cond, 2093 seebeck, 2094 kappa, 2095 hall, 2096 pn_doping_levels, 2097 mu_doping, 2098 seebeck_doping, 2099 cond_doping, 2100 kappa_doping, 2101 hall_doping, 2102 carrier_conc, 2103 ) = BoltztrapAnalyzer.parse_cond_and_hall(path_dir, doping_levels) 2104 2105 return BoltztrapAnalyzer( 2106 gap, 2107 mu_steps, 2108 cond, 2109 seebeck, 2110 kappa, 2111 hall, 2112 pn_doping_levels, 2113 mu_doping, 2114 seebeck_doping, 2115 cond_doping, 2116 kappa_doping, 2117 hall_doping, 2118 intrans, 2119 dos, 2120 pdos, 2121 carrier_conc, 2122 vol, 2123 warning, 2124 ) 2125 2126 if run_type == "DOS": 2127 trim = intrans["dos_type"] == "HISTO" 2128 dos, pdos = BoltztrapAnalyzer.parse_transdos(path_dir, efermi, dos_spin=dos_spin, trim_dos=trim) 2129 2130 return BoltztrapAnalyzer(gap=gap, dos=dos, dos_partial=pdos, warning=warning, vol=vol) 2131 2132 if run_type == "BANDS": 2133 bz_kpoints = np.loadtxt(os.path.join(path_dir, "boltztrap_band.dat"))[:, -3:] 2134 bz_bands = np.loadtxt(os.path.join(path_dir, "boltztrap_band.dat"))[:, 1:-6] 2135 return BoltztrapAnalyzer(bz_bands=bz_bands, bz_kpoints=bz_kpoints, warning=warning, vol=vol) 2136 2137 if run_type == "FERMI": 2138 """ """ 2139 2140 if os.path.exists(os.path.join(path_dir, "boltztrap_BZ.cube")): 2141 fs_data = read_cube_file(os.path.join(path_dir, "boltztrap_BZ.cube")) 2142 elif os.path.exists(os.path.join(path_dir, "fort.30")): 2143 fs_data = read_cube_file(os.path.join(path_dir, "fort.30")) 2144 else: 2145 raise BoltztrapError("No data file found for fermi surface") 2146 return BoltztrapAnalyzer(fermi_surface_data=fs_data) 2147 2148 raise ValueError("Run type: {} not recognized!".format(run_type)) 2149 2150 def as_dict(self): 2151 """ 2152 :return: MSONable dict. 2153 """ 2154 results = { 2155 "gap": self.gap, 2156 "mu_steps": self.mu_steps, 2157 "intrans": self.intrans, 2158 "cond": self._cond, 2159 "seebeck": self._seebeck, 2160 "kappa": self._kappa, 2161 "hall": self._hall, 2162 "doping": self.doping, 2163 "mu_doping": self.mu_doping, 2164 "seebeck_doping": self._seebeck_doping, 2165 "cond_doping": self._cond_doping, 2166 "kappa_doping": self._kappa_doping, 2167 "hall_doping": self._hall_doping, 2168 "dos": self.dos.as_dict(), 2169 "dos_partial": self._dos_partial, 2170 "carrier_conc": self._carrier_conc, 2171 "vol": self.vol, 2172 "warning": self.warning, 2173 } 2174 return jsanitize(results) 2175 2176 @staticmethod 2177 def from_dict(data): 2178 """ 2179 :param data: Dict representation. 2180 :return: BoltztrapAnalyzer 2181 """ 2182 2183 def _make_float_array(a): 2184 res = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] 2185 for i in range(3): 2186 for j in range(3): 2187 res[i][j] = float(a[i][j]) 2188 return res 2189 2190 def _make_float_hall(a): 2191 return list(a[:27]) 2192 2193 gap = data.get("gap") 2194 mu_steps = [float(d) for d in data["mu_steps"]] if "mu_steps" in data else None 2195 cond = ( 2196 {int(d): [_make_float_array(v) for v in data["cond"][d]] for d in data["cond"]} if "cond" in data else None 2197 ) 2198 seebeck = ( 2199 {int(d): [_make_float_array(v) for v in data["seebeck"][d]] for d in data["seebeck"]} 2200 if "seebeck" in data 2201 else None 2202 ) 2203 kappa = ( 2204 {int(d): [_make_float_array(v) for v in data["kappa"][d]] for d in data["kappa"]} 2205 if "kappa" in data 2206 else None 2207 ) 2208 hall = ( 2209 {int(d): [_make_float_hall(v) for v in data["hall"][d]] for d in data["hall"]} if "hall" in data else None 2210 ) 2211 doping = ( 2212 { 2213 "p": [float(d) for d in data["doping"]["p"]], 2214 "n": [float(d) for d in data["doping"]["n"]], 2215 } 2216 if "doping" in data 2217 else None 2218 ) 2219 2220 mu_doping = ( 2221 { 2222 "p": {int(d): [float(v) for v in data["mu_doping"]["p"][d]] for d in data["mu_doping"]["p"]}, 2223 "n": {int(d): [float(v) for v in data["mu_doping"]["n"][d]] for d in data["mu_doping"]["n"]}, 2224 } 2225 if "mu_doping" in data 2226 else None 2227 ) 2228 2229 seebeck_doping = ( 2230 { 2231 "p": { 2232 int(d): [_make_float_array(v) for v in data["seebeck_doping"]["p"][d]] 2233 for d in data["seebeck_doping"]["p"] 2234 }, 2235 "n": { 2236 int(d): [_make_float_array(v) for v in data["seebeck_doping"]["n"][d]] 2237 for d in data["seebeck_doping"]["n"] 2238 }, 2239 } 2240 if "seebeck_doping" in data 2241 else None 2242 ) 2243 2244 cond_doping = ( 2245 { 2246 "p": { 2247 int(d): [_make_float_array(v) for v in data["cond_doping"]["p"][d]] 2248 for d in data["cond_doping"]["p"] 2249 }, 2250 "n": { 2251 int(d): [_make_float_array(v) for v in data["cond_doping"]["n"][d]] 2252 for d in data["cond_doping"]["n"] 2253 }, 2254 } 2255 if "cond_doping" in data 2256 else None 2257 ) 2258 2259 kappa_doping = ( 2260 { 2261 "p": { 2262 int(d): [_make_float_array(v) for v in data["kappa_doping"]["p"][d]] 2263 for d in data["kappa_doping"]["p"] 2264 }, 2265 "n": { 2266 int(d): [_make_float_array(v) for v in data["kappa_doping"]["n"][d]] 2267 for d in data["kappa_doping"]["n"] 2268 }, 2269 } 2270 if "kappa_doping" in data 2271 else None 2272 ) 2273 2274 hall_doping = ( 2275 { 2276 "p": { 2277 int(d): [_make_float_hall(v) for v in data["hall_doping"]["p"][d]] for d in data["hall_doping"]["p"] 2278 }, 2279 "n": { 2280 int(d): [_make_float_hall(v) for v in data["hall_doping"]["n"][d]] for d in data["hall_doping"]["n"] 2281 }, 2282 } 2283 if "hall_doping" in data 2284 else None 2285 ) 2286 2287 dos = Dos.from_dict(data["dos"]) if "dos" in data else None 2288 dos_partial = data.get("dos_partial") 2289 carrier_conc = data.get("carrier_conc") 2290 vol = data.get("vol") 2291 warning = data.get("warning") 2292 2293 return BoltztrapAnalyzer( 2294 gap=gap, 2295 mu_steps=mu_steps, 2296 cond=cond, 2297 seebeck=seebeck, 2298 kappa=kappa, 2299 hall=hall, 2300 doping=doping, 2301 mu_doping=mu_doping, 2302 seebeck_doping=seebeck_doping, 2303 cond_doping=cond_doping, 2304 kappa_doping=kappa_doping, 2305 hall_doping=hall_doping, 2306 dos=dos, 2307 dos_partial=dos_partial, 2308 carrier_conc=carrier_conc, 2309 vol=vol, 2310 warning=warning, 2311 ) 2312 2313 2314def read_cube_file(filename): 2315 """ 2316 :param filename: Cube filename 2317 :return: Energy data. 2318 """ 2319 with open(filename, "rt") as f: 2320 natoms = 0 2321 count_line = 0 2322 for line in f: 2323 line = line.rstrip("\n") 2324 if count_line == 0 and "CUBE" not in line: 2325 raise ValueError("CUBE file format not recognized") 2326 2327 if count_line == 2: 2328 tokens = line.split() 2329 natoms = int(tokens[0]) 2330 if count_line == 3: 2331 tokens = line.split() 2332 n1 = int(tokens[0]) 2333 elif count_line == 4: 2334 tokens = line.split() 2335 n2 = int(tokens[0]) 2336 elif count_line == 5: 2337 tokens = line.split() 2338 n3 = int(tokens[0]) 2339 elif count_line > 5: 2340 break 2341 2342 count_line += 1 2343 2344 if "fort.30" in filename: 2345 energy_data = np.genfromtxt(filename, skip_header=natoms + 6, skip_footer=1) 2346 nlines_data = len(energy_data) 2347 last_line = np.genfromtxt(filename, skip_header=nlines_data + natoms + 6) 2348 energy_data = np.append(energy_data.flatten(), last_line).reshape(n1, n2, n3) 2349 elif "boltztrap_BZ.cube" in filename: 2350 energy_data = np.loadtxt(filename, skiprows=natoms + 6).reshape(n1, n2, n3) 2351 2352 energy_data /= Energy(1, "eV").to("Ry") 2353 2354 return energy_data 2355 2356 2357def compare_sym_bands(bands_obj, bands_ref_obj, nb=None): 2358 """ 2359 Compute the mean of correlation between bzt and vasp bandstructure on 2360 sym line, for all bands and locally (for each branches) the difference 2361 squared (%) if nb is specified. 2362 """ 2363 2364 if bands_ref_obj.is_spin_polarized: 2365 nbands = min(bands_obj.nb_bands, 2 * bands_ref_obj.nb_bands) 2366 else: 2367 # TODO: why is this needed? Shouldn't pmg take care of nb_bands? 2368 nbands = min(len(bands_obj.bands[Spin.up]), len(bands_ref_obj.bands[Spin.up])) 2369 # print(nbands) 2370 arr_bands = np.array(bands_obj.bands[Spin.up][:nbands]) 2371 # arr_bands_lavg = (arr_bands-np.mean(arr_bands,axis=1).reshape(nbands,1)) 2372 2373 if bands_ref_obj.is_spin_polarized: 2374 arr_bands_ref_up = np.array(bands_ref_obj.bands[Spin.up]) 2375 arr_bands_ref_dw = np.array(bands_ref_obj.bands[Spin.down]) 2376 # print(arr_bands_ref_up.shape) 2377 arr_bands_ref = np.vstack((arr_bands_ref_up, arr_bands_ref_dw)) 2378 arr_bands_ref = np.sort(arr_bands_ref, axis=0)[:nbands] 2379 # print(arr_bands_ref.shape) 2380 else: 2381 arr_bands_ref = np.array(bands_ref_obj.bands[Spin.up][:nbands]) 2382 2383 # arr_bands_ref_lavg = 2384 # (arr_bands_ref-np.mean(arr_bands_ref,axis=1).reshape(nbands,1)) 2385 2386 # err = np.sum((arr_bands_lavg-arr_bands_ref_lavg)**2,axis=1)/nkpt 2387 corr = np.array([distance.correlation(arr_bands[idx], arr_bands_ref[idx]) for idx in range(nbands)]) 2388 2389 if isinstance(nb, int): 2390 nb = [nb] 2391 2392 bcheck = {} 2393 2394 if max(nb) < nbands: 2395 branches = [[s["start_index"], s["end_index"], s["name"]] for s in bands_ref_obj.branches] 2396 2397 if not bands_obj.is_metal() and not bands_ref_obj.is_metal(): 2398 zero_ref = bands_ref_obj.get_vbm()["energy"] 2399 zero = bands_obj.get_vbm()["energy"] 2400 if not zero: 2401 vbm = bands_ref_obj.get_vbm()["band_index"][Spin.up][-1] 2402 zero = max(arr_bands[vbm]) 2403 else: 2404 zero_ref = 0 # bands_ref_obj.efermi 2405 zero = 0 # bands_obj.efermi 2406 print(zero, zero_ref) 2407 2408 for nbi in nb: 2409 bcheck[nbi] = {} 2410 2411 bcheck[nbi]["Dist"] = np.mean(abs(arr_bands[nbi] - zero - arr_bands_ref[nbi] + zero_ref)) 2412 bcheck[nbi]["Corr"] = corr[nbi] 2413 2414 for start, end, name in branches: 2415 # werr.append((sum((arr_bands_corr[nb][start:end+1] - 2416 # arr_bands_ref_corr[nb][start:end+1])**2)/(end+1-start)*100,name)) 2417 bcheck[nbi][name] = np.mean( 2418 abs(arr_bands[nbi][start : end + 1] - zero - arr_bands_ref[nbi][start : end + 1] + zero_ref) 2419 ) 2420 else: 2421 bcheck = "No nb given" 2422 2423 return bcheck 2424 2425 2426def seebeck_spb(eta, Lambda=0.5): 2427 """ 2428 Seebeck analytic formula in the single parabolic model 2429 """ 2430 from fdint import fdk 2431 2432 return ( 2433 constants.k 2434 / constants.e 2435 * ((2.0 + Lambda) * fdk(1.0 + Lambda, eta) / ((1.0 + Lambda) * fdk(Lambda, eta)) - eta) 2436 * 1e6 2437 ) 2438 2439 2440def eta_from_seebeck(seeb, Lambda): 2441 """ 2442 It takes a value of seebeck and adjusts the analytic seebeck until it's equal 2443 Returns: eta where the two seebeck coefficients are equal 2444 (reduced chemical potential) 2445 """ 2446 from scipy.optimize import fsolve 2447 2448 out = fsolve(lambda x: (seebeck_spb(x, Lambda) - abs(seeb)) ** 2, 1.0, full_output=True) 2449 return out[0][0] 2450 2451 2452def seebeck_eff_mass_from_carr(eta, n, T, Lambda): 2453 """ 2454 Calculate seebeck effective mass at a certain carrier concentration 2455 eta in kB*T units, n in cm-3, T in K, returns mass in m0 units 2456 """ 2457 try: 2458 from fdint import fdk 2459 except ImportError: 2460 raise BoltztrapError( 2461 "fdint module not found. Please, install it.\n" + "It is needed to calculate Fermi integral quickly." 2462 ) 2463 2464 return (2 * np.pi ** 2 * abs(n) * 10 ** 6 / (fdk(0.5, eta))) ** (2.0 / 3) / ( 2465 2 * constants.m_e * constants.k * T / (constants.h / 2 / np.pi) ** 2 2466 ) 2467 2468 2469def seebeck_eff_mass_from_seebeck_carr(seeb, n, T, Lambda): 2470 """ 2471 Find the chemical potential where analytic and calculated seebeck are identical 2472 and then calculate the seebeck effective mass at that chemical potential and 2473 a certain carrier concentration n 2474 """ 2475 eta = eta_from_seebeck(seeb, Lambda) 2476 mass = seebeck_eff_mass_from_carr(eta, n, T, Lambda) 2477 return mass 2478