1# coding: utf-8 2# Copyright (c) Pymatgen Development Team. 3# Distributed under the terms of the MIT License. 4 5""" 6This module defines classes to represent crystal orbital Hamilton 7populations (COHP) and integrated COHP (ICOHP), but can also be used 8for crystal orbital overlap populations (COOP). 9""" 10 11import re 12import sys 13import warnings 14 15import numpy as np 16from monty.json import MSONable 17from scipy.interpolate import InterpolatedUnivariateSpline 18from pymatgen.core.sites import PeriodicSite 19from pymatgen.core.structure import Structure 20from pymatgen.electronic_structure.core import Orbital, Spin 21from pymatgen.io.lmto import LMTOCopl 22from pymatgen.io.lobster import Cohpcar 23from pymatgen.util.coord import get_linear_interpolated_value 24from pymatgen.util.num import round_to_sigfigs 25 26__author__ = "Marco Esters, Janine George" 27__copyright__ = "Copyright 2017, The Materials Project" 28__version__ = "0.2" 29__maintainer__ = "Marco Esters, Janine George" 30__email__ = "esters@uoregon.edu, janine.george@uclouvain.be" 31__date__ = "Dec 13, 2017" 32 33 34class Cohp(MSONable): 35 """ 36 Basic COHP object. 37 """ 38 39 def __init__(self, efermi, energies, cohp, are_coops=False, are_cobis=False, icohp=None): 40 """ 41 Args: 42 are_coops: Indicates whether this object describes COOPs. 43 are_cobis: Indicates whether this object describes COBIs. 44 efermi: Fermi energy. 45 energies: A sequence of energies. 46 cohp ({Spin: np.array}): representing the COHP for each spin. 47 icohp ({Spin: np.array}): representing the ICOHP for each spin. 48 """ 49 self.are_coops = are_coops 50 self.are_cobis = are_cobis 51 self.efermi = efermi 52 self.energies = np.array(energies) 53 self.cohp = cohp 54 self.icohp = icohp 55 56 def __repr__(self): 57 return self.__str__() 58 59 def __str__(self): 60 """ 61 Returns a string that can be easily plotted (e.g. using gnuplot). 62 """ 63 if self.are_coops: 64 cohpstring = "COOP" 65 elif self.are_cobis: 66 cohpstring = "COBI" 67 else: 68 cohpstring = "COHP" 69 header = ["Energy", cohpstring + "Up"] 70 data = [self.energies, self.cohp[Spin.up]] 71 if Spin.down in self.cohp: 72 header.append(cohpstring + "Down") 73 data.append(self.cohp[Spin.down]) 74 if self.icohp: 75 header.append("I" + cohpstring + "Up") 76 data.append(self.icohp[Spin.up]) 77 if Spin.down in self.cohp: 78 header.append("I" + cohpstring + "Down") 79 data.append(self.icohp[Spin.down]) 80 formatheader = "#" + " ".join(["{:15s}" for __ in header]) 81 formatdata = " ".join(["{:.5f}" for __ in header]) 82 stringarray = [formatheader.format(*header)] 83 for i, __ in enumerate(self.energies): 84 stringarray.append(formatdata.format(*[d[i] for d in data])) 85 return "\n".join(stringarray) 86 87 def as_dict(self): 88 """ 89 Json-serializable dict representation of COHP. 90 """ 91 d = { 92 "@module": self.__class__.__module__, 93 "@class": self.__class__.__name__, 94 "are_coops": self.are_coops, 95 "are_cobis": self.are_cobis, 96 "efermi": self.efermi, 97 "energies": self.energies.tolist(), 98 "COHP": {str(spin): pops.tolist() for spin, pops in self.cohp.items()}, 99 } 100 if self.icohp: 101 d["ICOHP"] = {str(spin): pops.tolist() for spin, pops in self.icohp.items()} 102 return d 103 104 def get_cohp(self, spin=None, integrated=False): 105 """ 106 Returns the COHP or ICOHP for a particular spin. 107 108 Args: 109 spin: Spin. Can be parsed as spin object, integer (-1/1) 110 or str ("up"/"down") 111 integrated: Return COHP (False) or ICOHP (True) 112 113 Returns: 114 Returns the CHOP or ICOHP for the input spin. If Spin is 115 None and both spins are present, both spins will be returned 116 as a dictionary. 117 """ 118 if not integrated: 119 populations = self.cohp 120 else: 121 populations = self.icohp 122 123 if populations is None: 124 return None 125 if spin is None: 126 return populations 127 if isinstance(spin, int): 128 spin = Spin(spin) 129 elif isinstance(spin, str): 130 s = {"up": 1, "down": -1}[spin.lower()] 131 spin = Spin(s) 132 return {spin: populations[spin]} 133 134 def get_icohp(self, spin=None): 135 """ 136 Convenient alternative to get the ICOHP for a particular spin. 137 """ 138 return self.get_cohp(spin=spin, integrated=True) 139 140 def get_interpolated_value(self, energy, integrated=False): 141 """ 142 Returns the COHP for a particular energy. 143 144 Args: 145 energy: Energy to return the COHP value for. 146 """ 147 inter = {} 148 for spin in self.cohp: 149 if not integrated: 150 inter[spin] = get_linear_interpolated_value(self.energies, self.cohp[spin], energy) 151 elif self.icohp is not None: 152 inter[spin] = get_linear_interpolated_value(self.energies, self.icohp[spin], energy) 153 else: 154 raise ValueError("ICOHP is empty.") 155 return inter 156 157 def has_antibnd_states_below_efermi(self, spin=None, limit=0.01): 158 """ 159 Returns dict indicating if there are antibonding states below the Fermi level depending on the spin 160 spin: Spin 161 limit: -COHP smaller -limit will be considered. 162 163 """ 164 warnings.warn("This method has not been tested on many examples. Check the parameter limit, pls!") 165 166 populations = self.cohp 167 number_energies_below_efermi = len([x for x in self.energies if x <= self.efermi]) 168 169 if populations is None: 170 return None 171 if spin is None: 172 dict_to_return = {} 173 for sp, cohpvalues in populations.items(): 174 if (max(cohpvalues[0:number_energies_below_efermi])) > limit: 175 dict_to_return[sp] = True 176 else: 177 dict_to_return[sp] = False 178 else: 179 dict_to_return = {} 180 if isinstance(spin, int): 181 spin = Spin(spin) 182 elif isinstance(spin, str): 183 s = {"up": 1, "down": -1}[spin.lower()] 184 spin = Spin(s) 185 if (max(populations[spin][0:number_energies_below_efermi])) > limit: 186 dict_to_return[spin] = True 187 else: 188 dict_to_return[spin] = False 189 190 return dict_to_return 191 192 @classmethod 193 def from_dict(cls, d): 194 """ 195 Returns a COHP object from a dict representation of the COHP. 196 197 """ 198 if "ICOHP" in d: 199 icohp = {Spin(int(key)): np.array(val) for key, val in d["ICOHP"].items()} 200 else: 201 icohp = None 202 if "are_cobis" not in d: 203 are_cobis = False 204 else: 205 are_cobis = d["are_cobis"] 206 return Cohp( 207 d["efermi"], 208 d["energies"], 209 {Spin(int(key)): np.array(val) for key, val in d["COHP"].items()}, 210 icohp=icohp, 211 are_coops=d["are_coops"], 212 are_cobis=are_cobis, 213 ) 214 215 216class CompleteCohp(Cohp): 217 """ 218 A wrapper class that defines an average COHP, and individual COHPs. 219 220 .. attribute: are_coops 221 222 Indicates whether the object is consisting of COOPs. 223 224 225 .. attribute: are_cobis 226 227 Indicates whether the object is consisting of COBIs. 228 229 230 .. attribute: efermi 231 232 Fermi energy 233 234 .. attribute: energies 235 236 Sequence of energies 237 238 .. attribute: structure 239 240 Structure associated with the COHPs. 241 242 .. attribute: cohp, icohp 243 244 The average COHP/ICOHP. 245 246 .. attribute: all_cohps 247 248 A dict of COHPs for individual bonds of the form {label: COHP} 249 250 .. attribute: orb_res_cohp 251 252 Orbital-resolved COHPs. 253 """ 254 255 def __init__( 256 self, 257 structure, 258 avg_cohp, 259 cohp_dict, 260 bonds=None, 261 are_coops=False, 262 are_cobis=False, 263 orb_res_cohp=None, 264 ): 265 """ 266 Args: 267 structure: Structure assosciated with this COHP. 268 avg_cohp: The average cohp as a COHP object. 269 cohps: A dict of COHP objects for individual bonds of the form 270 {label: COHP} 271 bonds: A dict containing information on the bonds of the form 272 {label: {key: val}}. The key-val pair can be any information 273 the user wants to put in, but typically contains the sites, 274 the bond length, and the number of bonds. If nothing is 275 supplied, it will default to an empty dict. 276 are_coops: indicates whether the Cohp objects are COOPs. 277 Defauls to False for COHPs. 278 are_cobis: indicates whether the Cohp objects are COBIs. 279 Defauls to False for COHPs. 280 orb_res_cohp: Orbital-resolved COHPs. 281 """ 282 if are_coops and are_cobis: 283 raise ValueError("You cannot have info about COOPs and COBIs in the same file.") 284 super().__init__( 285 avg_cohp.efermi, 286 avg_cohp.energies, 287 avg_cohp.cohp, 288 are_coops=are_coops, 289 are_cobis=are_cobis, 290 icohp=avg_cohp.icohp, 291 ) 292 self.structure = structure 293 self.are_coops = are_coops 294 self.are_cobis = are_cobis 295 self.all_cohps = cohp_dict 296 self.orb_res_cohp = orb_res_cohp 297 if bonds is None: 298 self.bonds = {label: {} for label in self.all_cohps.keys()} 299 else: 300 self.bonds = bonds 301 302 def __str__(self): 303 if self.are_coops: 304 return "Complete COOPs for " + str(self.structure) 305 if self.are_cobis: 306 return "Complete COBIs for " + str(self.structure) 307 return "Complete COHPs for " + str(self.structure) 308 309 def as_dict(self): 310 """ 311 Json-serializable dict representation of CompleteCohp. 312 """ 313 d = { 314 "@module": self.__class__.__module__, 315 "@class": self.__class__.__name__, 316 "are_coops": self.are_coops, 317 "are_cobis": self.are_cobis, 318 "efermi": self.efermi, 319 "structure": self.structure.as_dict(), 320 "energies": self.energies.tolist(), 321 "COHP": {"average": {str(spin): pops.tolist() for spin, pops in self.cohp.items()}}, 322 } 323 324 if self.icohp is not None: 325 d["ICOHP"] = {"average": {str(spin): pops.tolist() for spin, pops in self.icohp.items()}} 326 327 for label in self.all_cohps.keys(): 328 d["COHP"].update({label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].cohp.items()}}) 329 if self.all_cohps[label].icohp is not None: 330 if "ICOHP" not in d.keys(): 331 d["ICOHP"] = { 332 label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].icohp.items()} 333 } 334 else: 335 d["ICOHP"].update( 336 {label: {str(spin): pops.tolist() for spin, pops in self.all_cohps[label].icohp.items()}} 337 ) 338 if False in [bond_dict == {} for bond_dict in self.bonds.values()]: 339 d["bonds"] = { 340 bond: { 341 "length": self.bonds[bond]["length"], 342 "sites": [site.as_dict() for site in self.bonds[bond]["sites"]], 343 } 344 for bond in self.bonds 345 } 346 if self.orb_res_cohp: 347 orb_dict = {} 348 for label in self.orb_res_cohp: 349 orb_dict[label] = {} 350 for orbs in self.orb_res_cohp[label]: 351 cohp = {str(spin): pops.tolist() for spin, pops in self.orb_res_cohp[label][orbs]["COHP"].items()} 352 orb_dict[label][orbs] = {"COHP": cohp} 353 icohp = {str(spin): pops.tolist() for spin, pops in self.orb_res_cohp[label][orbs]["ICOHP"].items()} 354 orb_dict[label][orbs]["ICOHP"] = icohp 355 orbitals = [[orb[0], orb[1].name] for orb in self.orb_res_cohp[label][orbs]["orbitals"]] 356 orb_dict[label][orbs]["orbitals"] = orbitals 357 d["orb_res_cohp"] = orb_dict 358 359 return d 360 361 def get_cohp_by_label(self, label, summed_spin_channels=False): 362 """ 363 Get specific COHP object. 364 365 Args: 366 label: string (for newer Lobster versions: a number) 367 summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true 368 369 Returns: 370 Returns the COHP object to simplify plotting 371 """ 372 # TODO: possibly add option to only get one spin channel as Spin.down in the future! 373 if label.lower() == "average": 374 divided_cohp = self.cohp 375 divided_icohp = self.icohp 376 377 else: 378 divided_cohp = self.all_cohps[label].get_cohp(spin=None, integrated=False) 379 divided_icohp = self.all_cohps[label].get_icohp(spin=None) 380 381 if summed_spin_channels and Spin.down in self.cohp: 382 final_cohp = {} 383 final_icohp = {} 384 final_cohp[Spin.up] = np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0) 385 final_icohp[Spin.up] = np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) 386 else: 387 final_cohp = divided_cohp 388 final_icohp = divided_icohp 389 390 return Cohp( 391 efermi=self.efermi, 392 energies=self.energies, 393 cohp=final_cohp, 394 are_coops=self.are_coops, 395 are_cobis=self.are_cobis, 396 icohp=final_icohp, 397 ) 398 399 def get_summed_cohp_by_label_list(self, label_list, divisor=1, summed_spin_channels=False): 400 """ 401 Returns a COHP object that includes a summed COHP divided by divisor 402 403 Args: 404 label_list: list of labels for the COHP that should be included in the summed cohp 405 divisor: float/int, the summed cohp will be divided by this divisor 406 summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true 407 Returns: 408 Returns a COHP object including a summed COHP 409 """ 410 # check if cohps are spinpolarized or not 411 first_cohpobject = self.get_cohp_by_label(label_list[0]) 412 summed_cohp = first_cohpobject.cohp.copy() 413 summed_icohp = first_cohpobject.icohp.copy() 414 for label in label_list[1:]: 415 cohp_here = self.get_cohp_by_label(label) 416 summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp_here.cohp[Spin.up]], axis=0) 417 418 if Spin.down in summed_cohp: 419 summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp_here.cohp[Spin.down]], axis=0) 420 421 summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], cohp_here.icohp[Spin.up]], axis=0) 422 423 if Spin.down in summed_icohp: 424 summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp_here.icohp[Spin.down]], axis=0) 425 426 divided_cohp = {} 427 divided_icohp = {} 428 divided_cohp[Spin.up] = np.divide(summed_cohp[Spin.up], divisor) 429 divided_icohp[Spin.up] = np.divide(summed_icohp[Spin.up], divisor) 430 if Spin.down in summed_cohp: 431 divided_cohp[Spin.down] = np.divide(summed_cohp[Spin.down], divisor) 432 divided_icohp[Spin.down] = np.divide(summed_icohp[Spin.down], divisor) 433 434 if summed_spin_channels and Spin.down in summed_cohp: 435 final_cohp = {} 436 final_icohp = {} 437 final_cohp[Spin.up] = np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0) 438 final_icohp[Spin.up] = np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) 439 else: 440 final_cohp = divided_cohp 441 final_icohp = divided_icohp 442 443 return Cohp( 444 efermi=first_cohpobject.efermi, 445 energies=first_cohpobject.energies, 446 cohp=final_cohp, 447 are_coops=first_cohpobject.are_coops, 448 are_cobis=first_cohpobject.are_coops, 449 icohp=final_icohp, 450 ) 451 452 def get_summed_cohp_by_label_and_orbital_list( 453 self, label_list, orbital_list, divisor=1, summed_spin_channels=False 454 ): 455 """ 456 Returns a COHP object that includes a summed COHP divided by divisor 457 458 Args: 459 label_list: list of labels for the COHP that should be included in the summed cohp 460 orbital_list: list of orbitals for the COHPs that should be included in the summed cohp (same order as 461 label_list) 462 divisor: float/int, the summed cohp will be divided by this divisor 463 summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true 464 465 Returns: 466 Returns a COHP object including a summed COHP 467 """ 468 # check length of label_list and orbital_list: 469 if not len(label_list) == len(orbital_list): 470 raise ValueError("label_list and orbital_list don't have the same length!") 471 # check if cohps are spinpolarized or not 472 first_cohpobject = self.get_orbital_resolved_cohp(label_list[0], orbital_list[0]) 473 summed_cohp = first_cohpobject.cohp.copy() 474 summed_icohp = first_cohpobject.icohp.copy() 475 for ilabel, label in enumerate(label_list[1:], 1): 476 cohp_here = self.get_orbital_resolved_cohp(label, orbital_list[ilabel]) 477 summed_cohp[Spin.up] = np.sum([summed_cohp[Spin.up], cohp_here.cohp.copy()[Spin.up]], axis=0) 478 if Spin.down in summed_cohp: 479 summed_cohp[Spin.down] = np.sum([summed_cohp[Spin.down], cohp_here.cohp.copy()[Spin.down]], axis=0) 480 summed_icohp[Spin.up] = np.sum([summed_icohp[Spin.up], cohp_here.icohp.copy()[Spin.up]], axis=0) 481 if Spin.down in summed_icohp: 482 summed_icohp[Spin.down] = np.sum([summed_icohp[Spin.down], cohp_here.icohp.copy()[Spin.down]], axis=0) 483 484 divided_cohp = {} 485 divided_icohp = {} 486 divided_cohp[Spin.up] = np.divide(summed_cohp[Spin.up], divisor) 487 divided_icohp[Spin.up] = np.divide(summed_icohp[Spin.up], divisor) 488 if Spin.down in summed_cohp: 489 divided_cohp[Spin.down] = np.divide(summed_cohp[Spin.down], divisor) 490 divided_icohp[Spin.down] = np.divide(summed_icohp[Spin.down], divisor) 491 492 if summed_spin_channels and Spin.down in divided_cohp: 493 final_cohp = {} 494 final_icohp = {} 495 496 final_cohp[Spin.up] = np.sum([divided_cohp[Spin.up], divided_cohp[Spin.down]], axis=0) 497 final_icohp[Spin.up] = np.sum([divided_icohp[Spin.up], divided_icohp[Spin.down]], axis=0) 498 else: 499 final_cohp = divided_cohp 500 final_icohp = divided_icohp 501 502 return Cohp( 503 efermi=first_cohpobject.efermi, 504 energies=first_cohpobject.energies, 505 cohp=final_cohp, 506 are_coops=first_cohpobject.are_coops, 507 are_cobis=first_cohpobject.are_cobis, 508 icohp=final_icohp, 509 ) 510 511 def get_orbital_resolved_cohp(self, label, orbitals, summed_spin_channels=False): 512 """ 513 Get orbital-resolved COHP. 514 515 Args: 516 label: bond label (Lobster: labels as in ICOHPLIST/ICOOPLIST.lobster). 517 518 orbitals: The orbitals as a label, or list or tuple of the form 519 [(n1, orbital1), (n2, orbital2)]. Orbitals can either be str, 520 int, or Orbital. 521 522 summed_spin_channels: bool, will sum the spin channels and return the sum in Spin.up if true 523 524 Returns: 525 A Cohp object if CompleteCohp contains orbital-resolved cohp, 526 or None if it doesn't. 527 528 Note: It currently assumes that orbitals are str if they aren't the 529 other valid types. This is not ideal, but the easiest way to 530 avoid unicode issues between python 2 and python 3. 531 """ 532 if self.orb_res_cohp is None: 533 return None 534 if isinstance(orbitals, (list, tuple)): 535 cohp_orbs = [d["orbitals"] for d in self.orb_res_cohp[label].values()] 536 orbs = [] 537 for orbital in orbitals: 538 if isinstance(orbital[1], int): 539 orbs.append(tuple((orbital[0], Orbital(orbital[1])))) 540 elif isinstance(orbital[1], Orbital): 541 orbs.append(tuple((orbital[0], orbital[1]))) 542 elif isinstance(orbital[1], str): 543 orbs.append(tuple((orbital[0], Orbital[orbital[1]]))) 544 else: 545 raise TypeError("Orbital must be str, int, or Orbital.") 546 orb_index = cohp_orbs.index(orbs) 547 orb_label = list(self.orb_res_cohp[label].keys())[orb_index] 548 elif isinstance(orbitals, str): 549 orb_label = orbitals 550 else: 551 raise TypeError("Orbitals must be str, list, or tuple.") 552 try: 553 icohp = self.orb_res_cohp[label][orb_label]["ICOHP"] 554 except KeyError: 555 icohp = None 556 557 start_cohp = self.orb_res_cohp[label][orb_label]["COHP"] 558 start_icohp = icohp 559 560 if summed_spin_channels and Spin.down in start_cohp: 561 final_cohp = {} 562 final_icohp = {} 563 final_cohp[Spin.up] = np.sum([start_cohp[Spin.up], start_cohp[Spin.down]], axis=0) 564 if start_icohp is not None: 565 final_icohp[Spin.up] = np.sum([start_icohp[Spin.up], start_icohp[Spin.down]], axis=0) 566 else: 567 final_cohp = start_cohp 568 final_icohp = start_icohp 569 570 return Cohp( 571 self.efermi, 572 self.energies, 573 final_cohp, 574 icohp=final_icohp, 575 are_coops=self.are_coops, 576 are_cobis=self.are_cobis, 577 ) 578 579 @classmethod 580 def from_dict(cls, d): 581 """ 582 Returns CompleteCohp object from dict representation. 583 """ 584 cohp_dict = {} 585 efermi = d["efermi"] 586 energies = d["energies"] 587 structure = Structure.from_dict(d["structure"]) 588 if "bonds" in d.keys(): 589 bonds = { 590 bond: { 591 "length": d["bonds"][bond]["length"], 592 "sites": tuple(PeriodicSite.from_dict(site) for site in d["bonds"][bond]["sites"]), 593 } 594 for bond in d["bonds"] 595 } 596 else: 597 bonds = None 598 for label in d["COHP"]: 599 cohp = {Spin(int(spin)): np.array(d["COHP"][label][spin]) for spin in d["COHP"][label]} 600 try: 601 icohp = {Spin(int(spin)): np.array(d["ICOHP"][label][spin]) for spin in d["ICOHP"][label]} 602 except KeyError: 603 icohp = None 604 if label == "average": 605 avg_cohp = Cohp(efermi, energies, cohp, icohp=icohp) 606 else: 607 cohp_dict[label] = Cohp(efermi, energies, cohp, icohp=icohp) 608 609 if "orb_res_cohp" in d.keys(): 610 orb_cohp = {} 611 for label in d["orb_res_cohp"]: 612 orb_cohp[label] = {} 613 for orb in d["orb_res_cohp"][label]: 614 cohp = { 615 Spin(int(s)): np.array(d["orb_res_cohp"][label][orb]["COHP"][s], dtype=float) 616 for s in d["orb_res_cohp"][label][orb]["COHP"] 617 } 618 try: 619 icohp = { 620 Spin(int(s)): np.array(d["orb_res_cohp"][label][orb]["ICOHP"][s], dtype=float) 621 for s in d["orb_res_cohp"][label][orb]["ICOHP"] 622 } 623 except KeyError: 624 icohp = None 625 orbitals = [tuple((int(o[0]), Orbital[o[1]])) for o in d["orb_res_cohp"][label][orb]["orbitals"]] 626 orb_cohp[label][orb] = { 627 "COHP": cohp, 628 "ICOHP": icohp, 629 "orbitals": orbitals, 630 } 631 # If no total COHPs are present, calculate the total 632 # COHPs from the single-orbital populations. Total COHPs 633 # may not be present when the cohpgenerator keyword is used 634 # in LOBSTER versions 2.2.0 and earlier. 635 if label not in d["COHP"] or d["COHP"][label] is None: 636 cohp = { 637 Spin.up: np.sum( 638 np.array([orb_cohp[label][orb]["COHP"][Spin.up] for orb in orb_cohp[label]]), 639 axis=0, 640 ) 641 } 642 try: 643 cohp[Spin.down] = np.sum( 644 np.array([orb_cohp[label][orb]["COHP"][Spin.down] for orb in orb_cohp[label]]), 645 axis=0, 646 ) 647 except KeyError: 648 pass 649 650 orb_res_icohp = None in [orb_cohp[label][orb]["ICOHP"] for orb in orb_cohp[label]] 651 if (label not in d["ICOHP"] or d["ICOHP"][label] is None) and orb_res_icohp: 652 icohp = { 653 Spin.up: np.sum( 654 np.array([orb_cohp[label][orb]["ICOHP"][Spin.up] for orb in orb_cohp[label]]), 655 axis=0, 656 ) 657 } 658 try: 659 icohp[Spin.down] = np.sum( 660 np.array([orb_cohp[label][orb]["ICOHP"][Spin.down] for orb in orb_cohp[label]]), 661 axis=0, 662 ) 663 except KeyError: 664 pass 665 else: 666 orb_cohp = None 667 668 if "average" not in d["COHP"].keys(): 669 # calculate average 670 cohp = np.array([np.array(c) for c in d["COHP"].values()]).mean(axis=0) 671 try: 672 icohp = np.array([np.array(c) for c in d["ICOHP"].values()]).mean(axis=0) 673 except KeyError: 674 icohp = None 675 avg_cohp = Cohp(efermi, energies, cohp, icohp=icohp) 676 677 if "are_cobis" not in d: 678 are_cobis = False 679 else: 680 are_cobis = d["are_cobis"] 681 682 return CompleteCohp( 683 structure, 684 avg_cohp, 685 cohp_dict, 686 bonds=bonds, 687 are_coops=d["are_coops"], 688 are_cobis=are_cobis, 689 orb_res_cohp=orb_cohp, 690 ) 691 692 @classmethod 693 def from_file(cls, fmt, filename=None, structure_file=None, are_coops=False, are_cobis=False): 694 """ 695 Creates a CompleteCohp object from an output file of a COHP 696 calculation. Valid formats are either LMTO (for the Stuttgart 697 LMTO-ASA code) or LOBSTER (for the LOBSTER code). 698 699 Args: 700 cohp_file: Name of the COHP output file. Defaults to COPL 701 for LMTO and COHPCAR.lobster/COOPCAR.lobster for LOBSTER. 702 703 are_coops: Indicates whether the populations are COOPs or 704 COHPs. Defaults to False for COHPs. 705 706 are_cobis: Indicates whether the populations are COBIs or 707 COHPs. Defaults to False for COHPs. 708 709 fmt: A string for the code that was used to calculate 710 the COHPs so that the output file can be handled 711 correctly. Can take the values "LMTO" or "LOBSTER". 712 713 structure_file: Name of the file containing the structure. 714 If no file name is given, use CTRL for LMTO and POSCAR 715 for LOBSTER. 716 717 Returns: 718 A CompleteCohp object. 719 """ 720 if are_coops and are_cobis: 721 raise ValueError("You cannot have info about COOPs and COBIs in the same file.") 722 fmt = fmt.upper() 723 if fmt == "LMTO": 724 # LMTO COOPs and orbital-resolved COHP cannot be handled yet. 725 are_coops = False 726 are_cobis = False 727 orb_res_cohp = None 728 if structure_file is None: 729 structure_file = "CTRL" 730 if filename is None: 731 filename = "COPL" 732 cohp_file = LMTOCopl(filename=filename, to_eV=True) 733 elif fmt == "LOBSTER": 734 if are_coops and are_cobis: 735 raise ValueError("You cannot have info about COOPs and COBIs in the same file.") 736 if structure_file is None: 737 structure_file = "POSCAR" 738 if filename is None: 739 if filename is None: 740 if are_coops: 741 filename = "COOPCAR.lobster" 742 elif are_cobis: 743 filename = "COBICAR.lobster" 744 else: 745 filename = "COHPCAR.lobster" 746 cohp_file = Cohpcar(filename=filename, are_coops=are_coops, are_cobis=are_cobis) 747 orb_res_cohp = cohp_file.orb_res_cohp 748 else: 749 raise ValueError("Unknown format %s. Valid formats are LMTO " "and LOBSTER." % fmt) 750 751 structure = Structure.from_file(structure_file) 752 efermi = cohp_file.efermi 753 cohp_data = cohp_file.cohp_data 754 energies = cohp_file.energies 755 756 # Lobster shifts the energies so that the Fermi energy is at zero. 757 # Shifting should be done by the plotter object though. 758 759 spins = [Spin.up, Spin.down] if cohp_file.is_spin_polarized else [Spin.up] 760 if fmt == "LOBSTER": 761 energies += efermi 762 763 if orb_res_cohp is not None: 764 # If no total COHPs are present, calculate the total 765 # COHPs from the single-orbital populations. Total COHPs 766 # may not be present when the cohpgenerator keyword is used 767 # in LOBSTER versions 2.2.0 and earlier. 768 # TODO: Test this more extensively 769 # pylint: disable=E1133,E1136 770 for label in orb_res_cohp: 771 if cohp_file.cohp_data[label]["COHP"] is None: 772 # print(label) 773 cohp_data[label]["COHP"] = { 774 sp: np.sum( 775 [orb_res_cohp[label][orbs]["COHP"][sp] for orbs in orb_res_cohp[label]], 776 axis=0, 777 ) 778 for sp in spins 779 } 780 if cohp_file.cohp_data[label]["ICOHP"] is None: 781 cohp_data[label]["ICOHP"] = { 782 sp: np.sum( 783 [orb_res_cohp[label][orbs]["ICOHP"][sp] for orbs in orb_res_cohp[label]], 784 axis=0, 785 ) 786 for sp in spins 787 } 788 789 if fmt == "LMTO": 790 # Calculate the average COHP for the LMTO file to be 791 # consistent with LOBSTER output. 792 avg_data = {"COHP": {}, "ICOHP": {}} 793 for i in avg_data: # pylint: disable=C0206 794 for spin in spins: 795 rows = np.array([v[i][spin] for v in cohp_data.values()]) 796 avg = np.average(rows, axis=0) 797 # LMTO COHPs have 5 significant figures 798 avg_data[i].update({spin: np.array([round_to_sigfigs(a, 5) for a in avg], dtype=float)}) 799 avg_cohp = Cohp(efermi, energies, avg_data["COHP"], icohp=avg_data["ICOHP"]) 800 else: 801 avg_cohp = Cohp( 802 efermi, 803 energies, 804 cohp_data["average"]["COHP"], 805 icohp=cohp_data["average"]["COHP"], 806 are_coops=are_coops, 807 are_cobis=are_cobis, 808 ) 809 del cohp_data["average"] 810 811 cohp_dict = { 812 label: Cohp(efermi, energies, v["COHP"], icohp=v["ICOHP"], are_coops=are_coops, are_cobis=are_cobis) 813 for label, v in cohp_data.items() 814 } 815 816 bond_dict = { 817 label: { 818 "length": v["length"], 819 "sites": [structure.sites[site] for site in v["sites"]], 820 } 821 for label, v in cohp_data.items() 822 } 823 824 return CompleteCohp( 825 structure, 826 avg_cohp, 827 cohp_dict, 828 bonds=bond_dict, 829 are_coops=are_coops, 830 are_cobis=are_cobis, 831 orb_res_cohp=orb_res_cohp, 832 ) 833 834 835class IcohpValue(MSONable): 836 """ 837 Class to store information on an ICOHP or ICOOP value 838 839 .. attribute:: num_bonds 840 number of bonds used for the average cohp (relevant for Lobster versions <3.0) (int) 841 842 .. attribute:: are_coops 843 Boolean to indicates ICOOPs 844 845 .. attribute:: are_cobis 846 Boolean to indicates ICOBIs 847 848 .. attribute:: icohp 849 dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down} 850 851 .. attribute:: summed_icohp: 852 sum of icohp/icoop of both spin channels 853 854 """ 855 856 def __init__(self, label, atom1, atom2, length, translation, num, icohp, are_coops=False, are_cobis=False): 857 """ 858 Args: 859 label: label for the icohp 860 atom1: str of atom that is contributing to the bond 861 atom2: str of second atom that is contributing to the bond 862 length: float of bond lengths 863 translation: translation list, e.g. [0,0,0] 864 num: integer describing how often the bond exists 865 icohp: dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down} 866 are_coops: if True, this are COOPs 867 are_cobis: if True, this are COBIs 868 """ 869 if are_coops and are_cobis: 870 raise ValueError("You cannot have info about COOPs and COBIs in the same file.") 871 self._are_coops = are_coops 872 self._are_cobis = are_cobis 873 self._label = label 874 self._atom1 = atom1 875 self._atom2 = atom2 876 self._length = length 877 self._translation = translation 878 self._num = num 879 self._icohp = icohp 880 if Spin.down in self._icohp: 881 self._is_spin_polarized = True 882 else: 883 self._is_spin_polarized = False 884 885 def __str__(self): 886 887 if not self._are_coops and not self._are_cobis: 888 if self._is_spin_polarized: 889 return ( 890 "ICOHP " 891 + str(self._label) 892 + " between " 893 + str(self._atom1) 894 + " and " 895 + str(self._atom2) 896 + " (" 897 + str(self._translation) 898 + "): " 899 + str(self._icohp[Spin.up]) 900 + " eV (Spin up) and " 901 + str(self._icohp[Spin.down]) 902 + " eV (Spin down)" 903 ) 904 return ( 905 "ICOHP " 906 + str(self._label) 907 + " between " 908 + str(self._atom1) 909 + " and " 910 + str(self._atom2) 911 + " (" 912 + str(self._translation) 913 + "): " 914 + str(self._icohp[Spin.up]) 915 + " eV (Spin up)" 916 ) 917 if self._are_coops: 918 if self._is_spin_polarized: 919 return ( 920 "ICOOP " 921 + str(self._label) 922 + " between " 923 + str(self._atom1) 924 + " and " 925 + str(self._atom2) 926 + " (" 927 + str(self._translation) 928 + "): " 929 + str(self._icohp[Spin.up]) 930 + " (Spin up) and " 931 + str(self._icohp[Spin.down]) 932 + " (Spin down)" 933 ) 934 return ( 935 "ICOOP " 936 + str(self._label) 937 + " between " 938 + str(self._atom1) 939 + " and " 940 + str(self._atom2) 941 + " (" 942 + str(self._translation) 943 + "): " 944 + str(self._icohp[Spin.up]) 945 + " (Spin up)" 946 ) 947 if self._is_spin_polarized: 948 return ( 949 "ICOBI " 950 + str(self._label) 951 + " between " 952 + str(self._atom1) 953 + " and " 954 + str(self._atom2) 955 + " (" 956 + str(self._translation) 957 + "): " 958 + str(self._icohp[Spin.up]) 959 + " (Spin up) and " 960 + str(self._icohp[Spin.down]) 961 + " (Spin down)" 962 ) 963 return ( 964 "ICOBI " 965 + str(self._label) 966 + " between " 967 + str(self._atom1) 968 + " and " 969 + str(self._atom2) 970 + " (" 971 + str(self._translation) 972 + "): " 973 + str(self._icohp[Spin.up]) 974 + " (Spin up)" 975 ) 976 977 @property 978 def num_bonds(self): 979 """ 980 tells the number of bonds for which the ICOHP value is an average 981 Returns: 982 Int 983 """ 984 return self._num 985 986 @property 987 def are_coops(self): 988 """ 989 tells if ICOOPs or not 990 Returns: 991 Boolean 992 """ 993 return self._are_coops 994 995 @property 996 def are_cobis(self): 997 """ 998 tells if ICOBIs or not 999 Returns: 1000 Boolean 1001 """ 1002 return self._are_cobis 1003 1004 @property 1005 def is_spin_polarized(self): 1006 """ 1007 tells if spin polarized calculation or not 1008 Returns: 1009 Boolean 1010 1011 """ 1012 return self._is_spin_polarized 1013 1014 def icohpvalue(self, spin=Spin.up): 1015 """ 1016 Args: 1017 spin: Spin.up or Spin.down 1018 Returns: 1019 icohpvalue (float) corresponding to chosen spin 1020 """ 1021 if not self.is_spin_polarized and spin == Spin.down: 1022 raise ValueError("The calculation was not performed with spin polarization") 1023 1024 return self._icohp[spin] 1025 1026 @property 1027 def icohp(self): 1028 """ 1029 dict with icohps for spinup and spindown 1030 Return: 1031 dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down} 1032 """ 1033 return self._icohp 1034 1035 @property 1036 def summed_icohp(self): 1037 """ 1038 Adds ICOHPs of both spin channels for spin polarized compounds 1039 Returns: 1040 icohp value in eV 1041 """ 1042 if self._is_spin_polarized: 1043 sum_icohp = self._icohp[Spin.down] + self._icohp[Spin.up] 1044 else: 1045 sum_icohp = self._icohp[Spin.up] 1046 return sum_icohp 1047 1048 1049class IcohpCollection(MSONable): 1050 """ 1051 Class to store IcohpValues 1052 1053 .. attribute:: are_coops 1054 Boolean to indicate if this are ICOOPs 1055 1056 .. attribute:: are_cobis 1057 Boolean to indicate if this are ICOOPs 1058 1059 .. attribute:: is_spin_polarized 1060 Boolean to indicate if the Lobster calculation was done spin polarized or not 1061 1062 """ 1063 1064 def __init__( 1065 self, 1066 list_labels, 1067 list_atom1, 1068 list_atom2, 1069 list_length, 1070 list_translation, 1071 list_num, 1072 list_icohp, 1073 is_spin_polarized, 1074 are_coops=False, 1075 are_cobis=False, 1076 ): 1077 """ 1078 Args: 1079 is_spin_polarized: Boolean to indicate if the Lobster calculation was done spin polarized or not Boolean to 1080 indicate if the Lobster calculation was done spin polarized or not 1081 are_coops: Boolean to indicate whether ICOOPs are stored 1082 are_cobis: Boolean to indicate whether ICOBIs are stored 1083 list_labels: list of labels for ICOHP/ICOOP values 1084 list_atom1: list of str of atomnames e.g. "O1" 1085 list_atom2: list of str of atomnames e.g. "O1" 1086 list_length: list of lengths of corresponding bonds in Angstrom 1087 list_translation: list of translation list, e.g. [0,0,0] 1088 list_num: list of equivalent bonds, usually 1 starting from Lobster 3.0.0 1089 list_icohp: list of dict={Spin.up: icohpvalue for spin.up, Spin.down: icohpvalue for spin.down} 1090 """ 1091 if are_coops and are_cobis: 1092 raise ValueError("You cannot have info about COOPs and COBIs in the same file.") 1093 self._are_coops = are_coops 1094 self._are_cobis = are_cobis 1095 self._icohplist = {} 1096 self._is_spin_polarized = is_spin_polarized 1097 self._list_labels = list_labels 1098 self._list_atom1 = list_atom1 1099 self._list_atom2 = list_atom2 1100 self._list_length = list_length 1101 self._list_translation = list_translation 1102 self._list_num = list_num 1103 self._list_icohp = list_icohp 1104 1105 for ilist, listel in enumerate(list_labels): 1106 self._icohplist[listel] = IcohpValue( 1107 label=listel, 1108 atom1=list_atom1[ilist], 1109 atom2=list_atom2[ilist], 1110 length=list_length[ilist], 1111 translation=list_translation[ilist], 1112 num=list_num[ilist], 1113 icohp=list_icohp[ilist], 1114 are_coops=are_coops, 1115 are_cobis=are_cobis, 1116 ) 1117 1118 def __str__(self): 1119 joinstr = [] 1120 for value in self._icohplist.values(): 1121 joinstr.append(str(value)) 1122 return "\n".join(joinstr) 1123 1124 def get_icohp_by_label(self, label, summed_spin_channels=True, spin=Spin.up): 1125 """ 1126 get an icohp value for a certain bond as indicated by the label (bond labels starting by "1" as in 1127 ICOHPLIST/ICOOPLIST) 1128 1129 Args: 1130 label: label in str format (usually the bond number in Icohplist.lobster/Icooplist.lobster 1131 summed_spin_channels: Boolean to indicate whether the ICOHPs/ICOOPs of both spin channels should be summed 1132 spin: if summed_spin_channels is equal to False, this spin indicates which spin channel should be returned 1133 1134 Returns: 1135 float describing ICOHP/ICOOP value 1136 """ 1137 1138 icohp_here = self._icohplist[label] 1139 if icohp_here._is_spin_polarized: 1140 if summed_spin_channels: 1141 return icohp_here.summed_icohp 1142 return icohp_here.icohpvalue(spin) 1143 return icohp_here.icohpvalue(spin) 1144 1145 def get_summed_icohp_by_label_list(self, label_list, divisor=1.0, summed_spin_channels=True, spin=Spin.up): 1146 """ 1147 get the sum of several ICOHP values that are indicated by a list of labels (labels of the bonds are the same as 1148 in ICOHPLIST/ICOOPLIST) 1149 1150 Args: 1151 label_list: list of labels of the ICOHPs/ICOOPs that should be summed 1152 divisor: is used to divide the sum 1153 summed_spin_channels: Boolean to indicate whether the ICOHPs/ICOOPs of both spin channels should be summed 1154 spin: if summed_spin_channels is equal to False, this spin indicates which spin channel should be returned 1155 1156 Returns: 1157 float that is a sum of all ICOHPs/ICOOPs as indicated with label_list 1158 """ 1159 sum_icohp = 0 1160 for label in label_list: 1161 icohp_here = self._icohplist[label] 1162 if icohp_here.num_bonds != 1: 1163 warnings.warn("One of the ICOHP values is an average over bonds. This is currently not considered.") 1164 # prints warning if num_bonds is not equal to 1 1165 if icohp_here._is_spin_polarized: 1166 if summed_spin_channels: 1167 sum_icohp = sum_icohp + icohp_here.summed_icohp 1168 else: 1169 sum_icohp = sum_icohp + icohp_here.icohpvalue(spin) 1170 else: 1171 sum_icohp = sum_icohp + icohp_here.icohpvalue(spin) 1172 return sum_icohp / divisor 1173 1174 def get_icohp_dict_by_bondlengths(self, minbondlength=0.0, maxbondlength=8.0): 1175 """ 1176 get a dict of IcohpValues corresponding to certain bond lengths 1177 Args: 1178 minbondlength: defines the minimum of the bond lengths of the bonds 1179 maxbondlength: defines the maximum of the bond lengths of the bonds 1180 Returns: 1181 dict of IcohpValues, the keys correspond to the values from the initial list_labels 1182 """ 1183 newicohp_dict = {} 1184 for value in self._icohplist.values(): 1185 if value._length >= minbondlength and value._length <= maxbondlength: 1186 newicohp_dict[value._label] = value 1187 return newicohp_dict 1188 1189 def get_icohp_dict_of_site( 1190 self, 1191 site, 1192 minsummedicohp=None, 1193 maxsummedicohp=None, 1194 minbondlength=0.0, 1195 maxbondlength=8.0, 1196 only_bonds_to=None, 1197 ): 1198 """ 1199 get a dict of IcohpValue for a certain site (indicated by integer) 1200 Args: 1201 site: integer describing the site of interest, order as in Icohplist.lobster/Icooplist.lobster, starts at 0 1202 minsummedicohp: float, minimal icohp/icoop of the bonds that are considered. It is the summed ICOHP value 1203 from both spin channels for spin polarized cases 1204 maxsummedicohp: float, maximal icohp/icoop of the bonds that are considered. It is the summed ICOHP value 1205 from both spin channels for spin polarized cases 1206 minbondlength: float, defines the minimum of the bond lengths of the bonds 1207 maxbondlength: float, defines the maximum of the bond lengths of the bonds 1208 only_bonds_to: list of strings describing the bonding partners that are allowed, e.g. ['O'] 1209 Returns: 1210 dict of IcohpValues, the keys correspond to the values from the initial list_labels 1211 """ 1212 1213 newicohp_dict = {} 1214 for key, value in self._icohplist.items(): 1215 atomnumber1 = int(re.split(r"(\d+)", value._atom1)[1]) - 1 1216 atomnumber2 = int(re.split(r"(\d+)", value._atom2)[1]) - 1 1217 if site in (atomnumber1, atomnumber2): 1218 # manipulate order of atoms so that searched one is always atom1 1219 if site == atomnumber2: 1220 save = value._atom1 1221 value._atom1 = value._atom2 1222 value._atom2 = save 1223 1224 if only_bonds_to is None: 1225 second_test = True 1226 else: 1227 second_test = re.split(r"(\d+)", value._atom2)[0] in only_bonds_to 1228 if value._length >= minbondlength and value._length <= maxbondlength and second_test: 1229 if minsummedicohp is not None: 1230 if value.summed_icohp >= minsummedicohp: 1231 if maxsummedicohp is not None: 1232 if value.summed_icohp <= maxsummedicohp: 1233 newicohp_dict[key] = value 1234 else: 1235 newicohp_dict[key] = value 1236 else: 1237 if maxsummedicohp is not None: 1238 if value.summed_icohp <= maxsummedicohp: 1239 newicohp_dict[key] = value 1240 else: 1241 newicohp_dict[key] = value 1242 1243 return newicohp_dict 1244 1245 def extremum_icohpvalue(self, summed_spin_channels=True, spin=Spin.up): 1246 """ 1247 get ICOHP/ICOOP of strongest bond 1248 Args: 1249 summed_spin_channels: Boolean to indicate whether the ICOHPs/ICOOPs of both spin channels should be summed 1250 1251 spin: if summed_spin_channels is equal to False, this spin indicates which spin channel should be returned 1252 Returns: 1253 lowest ICOHP/largest ICOOP value (i.e. ICOHP/ICOOP value of strongest bond) 1254 """ 1255 if self._are_coops or self._are_cobis: 1256 extremum = -sys.float_info.max 1257 else: 1258 extremum = sys.float_info.max 1259 1260 if not self._is_spin_polarized: 1261 if spin == Spin.down: 1262 warnings.warn("This spin channel does not exist. I am switching to Spin.up") 1263 spin = Spin.up 1264 1265 for value in self._icohplist.values(): 1266 if not value.is_spin_polarized or not summed_spin_channels: 1267 if not self._are_coops and not self._are_cobis: 1268 if value.icohpvalue(spin) < extremum: 1269 extremum = value.icohpvalue(spin) 1270 # print(extremum) 1271 else: 1272 if value.icohpvalue(spin) > extremum: 1273 extremum = value.icohpvalue(spin) 1274 # print(extremum) 1275 else: 1276 if not self._are_coops and not self._are_cobis: 1277 if value.summed_icohp < extremum: 1278 extremum = value.summed_icohp 1279 # print(extremum) 1280 else: 1281 if value.summed_icohp > extremum: 1282 extremum = value.summed_icohp 1283 # print(extremum) 1284 return extremum 1285 1286 @property 1287 def is_spin_polarized(self): 1288 """ 1289 :return: Whether it is spin polarized. 1290 """ 1291 return self._is_spin_polarized 1292 1293 @property 1294 def are_coops(self): 1295 """ 1296 :return: Whether this is a coop. 1297 """ 1298 return self._are_coops 1299 1300 @property 1301 def are_cobis(self): 1302 """ 1303 :return: Whether this a cobi. 1304 """ 1305 return self._are_cobis 1306 1307 1308def get_integrated_cohp_in_energy_range( 1309 cohp, label, orbital=None, energy_range=None, relative_E_Fermi=True, summed_spin_channels=True 1310): 1311 """ 1312 Method that can integrate completecohp objects which include data on integrated COHPs 1313 Args: 1314 cohp: CompleteCOHP object 1315 label: label of the COHP data 1316 orbital: If not None, a orbital resolved integrated COHP will be returned 1317 energy_range: if None, returns icohp value at Fermi level; 1318 if float, integrates from this float up to the Fermi level; 1319 if [float,float], will integrate in between 1320 relative_E_Fermi: if True, energy scale with E_Fermi at 0 eV is chosen 1321 summed_spin_channels: if True, Spin channels will be summed 1322 1323 Returns: 1324 float indicating the integrated COHP if summed_spin_channels==True, otherwise dict of the following form { 1325 Spin.up:float, Spin.down:float} 1326 """ 1327 summedicohp = {} 1328 if orbital is None: 1329 icohps = cohp.all_cohps[label].get_icohp(spin=None) 1330 if summed_spin_channels and Spin.down in icohps: 1331 summedicohp[Spin.up] = icohps[Spin.up] + icohps[Spin.down] 1332 else: 1333 summedicohp = icohps 1334 else: 1335 icohps = cohp.get_orbital_resolved_cohp(label=label, orbitals=orbital).icohp 1336 if summed_spin_channels and Spin.down in icohps: 1337 summedicohp[Spin.up] = icohps[Spin.up] + icohps[Spin.down] 1338 else: 1339 summedicohp = icohps 1340 1341 if energy_range is None: 1342 1343 energies_corrected = cohp.energies - cohp.efermi 1344 spl_spinup = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.up], ext=0) 1345 1346 if not summed_spin_channels and Spin.down in icohps: 1347 spl_spindown = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.down], ext=0) 1348 return {Spin.up: spl_spinup(0.0), Spin.down: spl_spindown(0.0)} 1349 if summed_spin_channels: 1350 return spl_spinup(0.0) 1351 1352 return {Spin.up: spl_spinup(0.0)} 1353 1354 # returns icohp value at the Fermi level! 1355 if isinstance(energy_range, float): 1356 if relative_E_Fermi: 1357 energies_corrected = cohp.energies - cohp.efermi 1358 spl_spinup = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.up], ext=0) 1359 1360 if not summed_spin_channels and Spin.down in icohps: 1361 spl_spindown = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.down], ext=0) 1362 return { 1363 Spin.up: spl_spinup(0) - spl_spinup(energy_range), 1364 Spin.down: spl_spindown(0) - spl_spindown(energy_range), 1365 } 1366 if summed_spin_channels: 1367 return spl_spinup(0) - spl_spinup(energy_range) 1368 return {Spin.up: spl_spinup(0) - spl_spinup(energy_range)} 1369 1370 energies_corrected = cohp.energies 1371 spl_spinup = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.up], ext=0) 1372 1373 if not summed_spin_channels and Spin.down in icohps: 1374 spl_spindown = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.down], ext=0) 1375 return { 1376 Spin.up: spl_spinup(cohp.efermi) - spl_spinup(energy_range), 1377 Spin.down: spl_spindown(cohp.efermi) - spl_spindown(energy_range), 1378 } 1379 if summed_spin_channels: 1380 return spl_spinup(cohp.efermi) - spl_spinup(energy_range) 1381 return {Spin.up: spl_spinup(cohp.efermi) - spl_spinup(energy_range)} 1382 1383 if relative_E_Fermi: 1384 energies_corrected = cohp.energies - cohp.efermi 1385 else: 1386 energies_corrected = cohp.energies 1387 1388 spl_spinup = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.up], ext=0) 1389 1390 if not summed_spin_channels and Spin.down in icohps: 1391 spl_spindown = InterpolatedUnivariateSpline(energies_corrected, summedicohp[Spin.down], ext=0) 1392 return { 1393 Spin.up: spl_spinup(energy_range[1]) - spl_spinup(energy_range[0]), 1394 Spin.down: spl_spindown(energy_range[1]) - spl_spindown(energy_range[0]), 1395 } 1396 if summed_spin_channels: 1397 return spl_spinup(energy_range[1]) - spl_spinup(energy_range[0]) 1398 1399 return {Spin.up: spl_spinup(energy_range[1]) - spl_spinup(energy_range[0])} 1400