+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""Script to analyze Lobster outputs from the command line."""
+from __future__ import annotations
+
+import argparse
+import json
+from math import log, sqrt
+from pathlib import Path
+
+import matplotlib.style
+from pymatgen.electronic_structure.cohp import CompleteCohp
+from pymatgen.io.lobster import Icohplist
+
+from lobsterpy.cohp.analyze import Analysis
+from lobsterpy.cohp.describe import Description
+from lobsterpy.plotting import (
+ IcohpDistancePlotter,
+ PlainCohpPlotter,
+ PlainDosPlotter,
+ get_style_list,
+)
+
+
+[docs]def main() -> None:
+ """Entry point for setup.py installer."""
+ args = get_parser().parse_args()
+ run(args)
+
+
+[docs]def get_parser() -> argparse.ArgumentParser:
+ """Construct argumentparser with subcommands and sections."""
+ parser = argparse.ArgumentParser(
+ description="Analyze and plot results from Lobster runs."
+ )
+
+ # Arguments common to all actions
+ input_parent = argparse.ArgumentParser(add_help=False)
+ input_file_group = input_parent.add_argument_group("Input files")
+ input_file_group.add_argument(
+ "--poscar",
+ "--POSCAR",
+ dest="poscar",
+ default="POSCAR",
+ type=Path,
+ help='path to POSCAR. Default is "POSCAR"',
+ )
+ input_file_group.add_argument(
+ "--charge",
+ default="CHARGE.lobster",
+ type=Path,
+ help='path to CHARGE.lobster. Default is "CHARGE.lobster"',
+ )
+ input_file_group.add_argument(
+ "--icohplist",
+ default="ICOHPLIST.lobster",
+ type=Path,
+ help='path to ICOHPLIST.lobster. Default is "ICOHPLIST.lobster"',
+ )
+ input_file_group.add_argument(
+ "--cohpcar",
+ default="COHPCAR.lobster",
+ type=Path,
+ help='path to COHPCAR.lobster. Default is "COHPCAR.lobster". This argument'
+ "will also be read when COBICARs or COOPCARs are plotted.",
+ )
+ input_file_group.add_argument(
+ "--incar",
+ default="INCAR",
+ type=Path,
+ help='path to INCAR. Default is "INCAR".',
+ )
+ input_file_group.add_argument(
+ "--potcar",
+ default="POTCAR",
+ type=Path,
+ help='path to POTCAR. Default is "POTCAR".',
+ )
+ input_file_group.add_argument(
+ "--potcar-symbols",
+ dest="potcarsymbols",
+ type=_potcar_symbols,
+ help="List of potcar symbols",
+ )
+
+ input_file_group.add_argument(
+ "--lobsterin",
+ default="lobsterin",
+ type=Path,
+ help='path to lobsterin. Default is "lobsterin".',
+ )
+
+ input_file_group.add_argument(
+ "--lobsterout",
+ default="lobsterout",
+ type=Path,
+ help='path to lobsterout. Default is "lobsterout".',
+ )
+
+ input_file_group.add_argument(
+ "--doscar",
+ default="DOSCAR.lobster",
+ type=Path,
+ help='path to DOSCAR.lobster. Default is "DOSCAR.lobster".',
+ )
+
+ input_file_group.add_argument(
+ "--vasprun",
+ default="vasprun.xml",
+ type=Path,
+ help='path to vasprun.xml. Default is "vasprun.xml".',
+ )
+
+ input_file_group.add_argument(
+ "--bandoverlaps",
+ default="bandOverlaps.lobster",
+ type=Path,
+ help='path to bandOverlaps.lobster. Default is "bandOverlaps.lobster".',
+ )
+
+ output_parent = argparse.ArgumentParser(add_help=False)
+ output_file_group = output_parent.add_argument_group("Output files")
+ output_file_group.add_argument(
+ "--incar-out",
+ "--incarout",
+ dest="incarout",
+ default="INCAR.lobsterpy",
+ type=Path,
+ help='path to INCAR that will be generated by lobsterpy. Default is "INCAR.lobsterpy"',
+ )
+ output_file_group.add_argument(
+ "--lobsterin-out",
+ "--lobsterinout",
+ dest="lobsterinout",
+ default="lobsterin.lobsterpy",
+ type=Path,
+ help='base for path to lobsterins that will be generated by lobsterpy. Default is "lobsterin.lobsterpy"',
+ )
+ output_file_group.add_argument(
+ "--overwrite",
+ "--overwrite-files",
+ dest="overwrite",
+ default=False,
+ action="store_true",
+ help="overwrites already created INCARs an lobsterins with the give name.",
+ )
+ # TODO: Add some output arguments: options to supply your own basis
+ output_file_group.add_argument(
+ "--userbasis",
+ "--user-basis",
+ default=None,
+ type=_element_basis,
+ nargs="+",
+ help="This setting will rely on a specific basis provided by the user "
+ "(e.g., --userbasis Cr.3d.3p.4s N.2s.2p). Default is None.",
+ )
+
+ plotting_parent = argparse.ArgumentParser(add_help=False)
+ plotting_group = plotting_parent.add_argument_group("Plotting")
+
+ broadening_group = plotting_group.add_mutually_exclusive_group()
+ broadening_group.add_argument_group("Broadening")
+ broadening_group.add_argument(
+ "--sigma",
+ type=float,
+ default=None,
+ help="Standard deviation of Gaussian broadening.",
+ )
+
+ broadening_group.add_argument(
+ "--fwhm",
+ type=float,
+ default=None,
+ help="Full-width-half-maximum of Gaussian broadening.",
+ )
+
+ plotting_group.add_argument(
+ "--integrated",
+ action="store_true",
+ help="Show integrated cohp/cobi/coop plots.",
+ )
+
+ plotting_group.add_argument(
+ "--hideplot",
+ "--hide-plot",
+ action="store_true",
+ help="Will hide plots generated by automaticplot/auto-plot/autoplot/plot args",
+ )
+
+ plotting_group.add_argument(
+ "--ylim",
+ dest="ylim",
+ nargs=2,
+ default=None,
+ type=float,
+ help="Energy range for plots",
+ )
+ plotting_group.add_argument(
+ "--xlim",
+ dest="xlim",
+ nargs=2,
+ default=None,
+ type=float,
+ help="COHP/COBI/COOP range for plots",
+ )
+
+ plotting_group.add_argument(
+ "--style",
+ type=str,
+ nargs="+",
+ default=None,
+ help="Matplotlib style sheet(s) for plot appearance",
+ )
+ plotting_group.add_argument(
+ "--no-base-style",
+ "--nobasestyle",
+ action="store_true",
+ dest="no_base_style",
+ help=(
+ "Disable inbuilt style entirely. This may prevent interference with external "
+ "stylesheets when using --style."
+ ),
+ )
+ plotting_group.add_argument("--title", type=str, default="", help="Plot title")
+ plotting_group.add_argument(
+ "--save-plot",
+ "--saveplot",
+ "-s",
+ type=Path,
+ metavar="FILENAME",
+ default=None,
+ dest="save_plot",
+ help="Save plot to file",
+ )
+ plotting_group.add_argument(
+ "--width", type=float, default=None, help="Plot width in inches"
+ )
+ plotting_group.add_argument(
+ "--height", type=float, default=None, help="Plot height in inches"
+ )
+ plotting_group.add_argument(
+ "--fontsize", "--font-size", type=float, default=None, help="Base font size"
+ )
+ group = plotting_group.add_mutually_exclusive_group()
+ group.add_argument(
+ "--spddos",
+ "--spd-dos",
+ action="store_true",
+ help="Will add spd projected dos to the DOS plot",
+ )
+ group.add_argument(
+ "--elementdos",
+ "--el-dos",
+ action="store_true",
+ help="Will add DOS projections for each element to the DOS plot",
+ )
+ plotting_group.add_argument(
+ "--summedspins",
+ "--summed-spins",
+ action="store_true",
+ help="Will plot summed DOS",
+ )
+ plotting_group.add_argument(
+ "--element",
+ "--el",
+ type=str,
+ nargs="+",
+ default=None,
+ help="Will add spd DOS projections for requested element to the DOS plot",
+ )
+ plotting_group.add_argument(
+ "--orbital",
+ "--orb",
+ type=str,
+ nargs="+",
+ default=None,
+ help="Orbital name for the site for which DOS are to be added",
+ )
+ plotting_group.add_argument(
+ "--site",
+ type=int,
+ nargs="+",
+ default=None,
+ help="Site index in the crystal structure for which DOS need to be added",
+ )
+
+ advanced_plotting_args = argparse.ArgumentParser(add_help=False)
+ advanced_plotting_args.add_argument(
+ "--invertaxis",
+ "--invert-axis",
+ action="store_true",
+ help="Will invert plot axis of DOS or COOPs COHPs or COBIS",
+ )
+ advanced_plotting_args.add_argument(
+ "--addtotaldos",
+ "--add-total-dos",
+ action="store_true",
+ help="Will all total dos to the DOS plot",
+ )
+
+ auto_parent = argparse.ArgumentParser(add_help=False)
+ auto_group = auto_parent.add_argument_group("Automatic analysis")
+ auto_group.add_argument(
+ "--json",
+ nargs="?",
+ type=Path,
+ default=None,
+ metavar="FILENAME",
+ const=Path("lobsterpy.json"),
+ help="Write a JSON file with the most important information",
+ )
+ auto_group.add_argument(
+ "--allbonds",
+ "--all-bonds",
+ action="store_true",
+ default=False,
+ help="This option will force the automatc analysis to consider"
+ " all bonds, not only cation-anion bonds (default) ",
+ )
+ auto_group.add_argument(
+ "--noisecutoff",
+ type=float,
+ default=None,
+ help="Sets the lower limit of icohps or icoops or icobis considered in"
+ " automatic analysis",
+ )
+ auto_group.add_argument(
+ "--cutofficohp",
+ type=float,
+ default=0.1,
+ help="Only bonds that are stronger than cutoff_icoxx *strongest ICOHP "
+ " (ICOBI or ICOOP) will be considered for automatic analysis.",
+ )
+ auto_group.add_argument(
+ "--orbitalresolved",
+ "--orbital-resolved",
+ action="store_true",
+ default=False,
+ help="Will switch on orbital resolved analysis of (I)COHPs or (I)COBIs or (I)COOPs with all relevant orbitals.",
+ )
+ auto_group.add_argument(
+ "--orbitalcutoff",
+ "--orbital-cutoff",
+ type=float,
+ default=0.05,
+ help="Will only work when orbital wise analysis is switched on (--orbitalresolved) "
+ "and only orbital interactions that are stronger than orbitalintcutoff * 100 of relevant "
+ "bonds (ICOHP or ICOBI or ICOOP) will be considered in automatic analysis.",
+ )
+ analysis_switch = argparse.ArgumentParser(add_help=False)
+ analysis_group = analysis_switch.add_argument_group(
+ "Switch type of analysis or plots"
+ )
+ analysis_group.add_argument(
+ "--cobis",
+ "--cobis",
+ action="store_true",
+ help="This option will start automatic bonding analysis of COBIs or"
+ " plot COBIs",
+ )
+ analysis_group.add_argument(
+ "--coops",
+ "--coops",
+ action="store_true",
+ help="This option will start automatic bonding analysis of COOPs or"
+ " plot COOPs",
+ )
+
+ interactive_plotter_args = argparse.ArgumentParser(add_help=False)
+ interactive_plotter_group = interactive_plotter_args.add_argument_group(
+ "Options specific to interactive plotter"
+ )
+
+ interactive_plotter_group.add_argument(
+ "--labelresolved",
+ "--label-resolved",
+ action="store_true",
+ help="Will create automatic interactive plots with all relevant bond labels. "
+ "If not set, plots will consists of summed cohps. (This argument works only"
+ "for interactive plots) ",
+ )
+ interactive_plotter_group.add_argument(
+ "--orbitalplot",
+ "--orbital-plot",
+ action="store_true",
+ help="Will generate automatic interactive (I)COHP or (I)COBI or (I)COOP plots with all relevant orbitals "
+ "If used along with --labelresolved arg, plots will be further label resolved else,"
+ "plots will consists of summed orbital cohps. ",
+ )
+
+ auto_group.add_argument(
+ "--calcqualityjson",
+ nargs="?",
+ type=Path,
+ default=None,
+ metavar="FILENAME",
+ const=Path("calc_quality_json.json"),
+ help="Write a JSON file with the LOBSTER calc quality analysis",
+ )
+
+ auto_group.add_argument(
+ "--erange",
+ dest="erange",
+ nargs=2,
+ default=[-5, 0],
+ type=int,
+ help="Energy range for DOS comparisons",
+ )
+
+ auto_group.add_argument(
+ "--nbins",
+ dest="nbins",
+ default=None,
+ type=int,
+ help="Number of bins for DOS comparisons",
+ )
+
+ auto_group.add_argument(
+ "--doscomp",
+ "--dos-comp",
+ action="store_true",
+ default=False,
+ help="This option will force the DOS comparison in automatic LOBSTER calc quality analysis ",
+ )
+
+ auto_group.add_argument(
+ "--bvacomp",
+ "--bva-comp",
+ action="store_true",
+ default=False,
+ help="This option will force the BVA charge comparison in automatic LOBSTER calc quality analysis ",
+ )
+
+ subparsers = parser.add_subparsers(
+ dest="action",
+ required=True,
+ help="Use -h/--help after the chosen subcommand to see further options.",
+ )
+ subparsers.add_parser(
+ "description",
+ parents=[input_parent, auto_parent, analysis_switch],
+ help=(
+ "Deliver a text description of the COHPs or COBIS or COOP results from Lobster "
+ "and VASP"
+ ),
+ )
+ subparsers.add_parser(
+ "calc-description",
+ parents=[input_parent, auto_parent],
+ help=(
+ "Deliver a text description of the LOBSTER calc quality analysis. "
+ "Mandatory required files: POSCAR, POTCAR, lobsterout, lobsterin. "
+ "Optional files (BVA comparison): CHARGE.lobster, "
+ "(DOS comparison): DOSCAR.lobster, Vasprun.xml."
+ ),
+ )
+
+ subparsers.add_parser(
+ "automatic-plot",
+ aliases=["automaticplot", "auto-plot", "autoplot"],
+ parents=[input_parent, auto_parent, plotting_parent, analysis_switch],
+ help=(
+ "Plot most important COHPs or COBIs or COOPs automatically."
+ " This option also includes an automatic description."
+ ),
+ )
+
+ subparsers.add_parser(
+ "automatic-plot-ia",
+ aliases=["automaticplotia", "autoplotia", "auto-plot-ia"],
+ parents=[
+ input_parent,
+ auto_parent,
+ plotting_parent,
+ interactive_plotter_args,
+ analysis_switch,
+ ],
+ help=(
+ "Creates an interactive plot of most important COHPs or COBIs or COOPs automatically."
+ ),
+ )
+
+ subparsers.add_parser(
+ "create-inputs",
+ aliases=["createinputs"],
+ parents=[input_parent, output_parent],
+ help=(
+ "Will create inputs for lobster computation. It only works with PBE POTCARs."
+ ),
+ )
+
+ subparsers.add_parser(
+ "plot-dos",
+ aliases=["plotdos"],
+ parents=[input_parent, plotting_parent, advanced_plotting_args],
+ help=("Will plot DOS from lobster computation."),
+ )
+ subparsers.add_parser(
+ "plot-icohps-distances",
+ aliases=["ploticohpsdistances"],
+ parents=[input_parent, plotting_parent, analysis_switch],
+ help=("Will plot icohps with respect to bond lengths"),
+ )
+
+ # Mode for normal plotting (without automatic detection of relevant COHPs)
+ plot_parser = subparsers.add_parser(
+ "plot",
+ parents=[input_parent, plotting_parent, analysis_switch],
+ help="Plot specific COHPs/COBIs/COOPs based on bond numbers.",
+ )
+
+ plot_parser.add_argument(
+ "bond_numbers",
+ nargs="+",
+ type=int,
+ help="List of bond numbers, determining COHPs/COBIs/COOPs to include in plot.",
+ )
+ plot_grouping = plot_parser.add_mutually_exclusive_group()
+ plot_grouping.add_argument(
+ "--summed",
+ action="store_true",
+ help="Show a summed COHP",
+ )
+ plot_grouping.add_argument(
+ "--orbitalwise",
+ dest="orbitalwise",
+ nargs="+",
+ default=None,
+ type=str,
+ help=(
+ "Plot cohps of specific orbitals. e.g. to plot 2s-2s interaction of "
+ 'bond with label 1, use "lobsterpy plot 1 --orbitalwise 2s-2s". '
+ 'To plot all orbitalwise cohps of one bond, you can use "all" instead of "2s-2s". '
+ "To plot orbitalwise interactions of more than one bond, use, for example, "
+ '"lobsterpy plot 1 1 --orbitalwise "3s-3s" "2px-3s"'
+ ),
+ )
+ return parser
+
+
+[docs]def _element_basis(string: str):
+ """
+ Parse element and basis from string.
+
+ Args:
+ string: string to parse
+
+ Returns:
+ element, basis
+ """
+ cut_list = string.split(".")
+ element = cut_list[0]
+ basis = " ".join(cut_list[1:])
+ return element, basis
+
+
+[docs]def _potcar_symbols(string: str):
+ """
+ Parse string of potcar symbols and return a list.
+
+ Args:
+ string: string of potcar symbols
+
+ Returns:
+ list of potcar symbols
+ """
+ return string.split(" ")
+
+
+[docs]def _user_figsize(width, height, aspect=None):
+ """Get figsize options from user input, if any.
+
+ If only width x or height is provided, use a target aspect ratio to derive
+ the other one.
+
+ Returns a dict which can be merged into style kwargs
+ """
+ if width is None and height is None:
+ return {}
+ if width is not None and height is not None:
+ return {"figure.figsize": (width, height)}
+
+ if aspect is None:
+ aspect = (sqrt(5) + 1) / 2 # Golden ratio
+ if width is None:
+ return {"figure.figsize": (height * aspect, height)}
+ return {"figure.figsize": (width, width / aspect)}
+
+
+# TODO: add automatic functionality for COBIs, COOPs
+[docs]def run(args):
+ """
+ Run actions based on args.
+
+ Args:
+ args: args for cli
+
+ """
+ if args.action in [
+ "automaticplot",
+ "autoplot",
+ "auto-plot",
+ "description",
+ "automatic-plot",
+ "plot",
+ "automatic-plot-ia",
+ "interactplot",
+ "automaticplotia",
+ "auto-plot-ia",
+ "autoplotia",
+ ]:
+ # Check for .gz files exist for default values and update accordingly
+ default_files = {
+ "poscar": "POSCAR",
+ "charge": "CHARGE.lobster",
+ "icohplist": "ICOHPLIST.lobster",
+ "cohpcar": "COHPCAR.lobster",
+ }
+
+ for arg_name in default_files:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "Files necessary for automatic analysis of LOBSTER outputs "
+ "not found in the current directory"
+ )
+
+ if args.coops:
+ if args.action != "plot":
+ default_files_coops = {
+ "icohplist": "ICOOPLIST.lobster",
+ "cohpcar": "COOPCAR.lobster",
+ }
+ else:
+ default_files_coops = {
+ "cohpcar": "COOPCAR.lobster",
+ }
+ for arg_name, file_name in default_files_coops.items():
+ setattr(args, arg_name, Path(file_name))
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "Files required for automatic analysis of COOPs (ICOOPLIST.lobster and"
+ " COOPCAR.lobster) not found in the directory"
+ )
+
+ if args.cobis:
+ if args.action != "plot":
+ default_files_cobis = {
+ "icohplist": "ICOBILIST.lobster",
+ "cohpcar": "COBICAR.lobster",
+ }
+ else:
+ default_files_cobis = {
+ "cohpcar": "COBICAR.lobster",
+ }
+ for arg_name, file_name in default_files_cobis.items():
+ setattr(args, arg_name, Path(file_name))
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "Files required for automatic analysis of COOPs (ICOBILIST.lobster and"
+ " COBICAR.lobster) not found in the directory"
+ )
+
+ if args.action in ["automaticplot", "autoplot", "auto-plot"]:
+ args.action = "automatic-plot"
+
+ if args.action in ["automaticplotia", "automatic-plot-ia", "auto-plot-ia"]:
+ args.action = "automaticplotia"
+
+ if args.action in [
+ "description",
+ "automatic-plot",
+ "automatic-plot-ia",
+ "automaticplotia",
+ "auto-plot-ia",
+ "autoplotia",
+ ]:
+ which_bonds = "all" if args.allbonds else "cation-anion"
+
+ analyse = Analysis(
+ path_to_poscar=args.poscar,
+ path_to_charge=args.charge,
+ path_to_cohpcar=args.cohpcar,
+ path_to_icohplist=args.icohplist,
+ which_bonds=which_bonds,
+ are_coops=args.coops,
+ are_cobis=args.cobis,
+ noise_cutoff=args.noisecutoff,
+ cutoff_icohp=args.cutofficohp,
+ orbital_cutoff=args.orbitalcutoff,
+ orbital_resolved=args.orbitalresolved,
+ )
+
+ describe = Description(analysis_object=analyse)
+ describe.write_description()
+
+ if args.json is not None:
+ analysedict = analyse.condensed_bonding_analysis
+ with open(args.json, "w") as fd:
+ json.dump(analysedict, fd)
+
+ if args.action in [
+ "plot",
+ "automatic-plot",
+ "automaticplotia",
+ "automatic-plot-ia",
+ "auto-plot-ia",
+ "autoplotia",
+ "plot-dos",
+ "plotdos",
+ "plot-icohps-distances",
+ "ploticohpsdistances",
+ ]:
+ style_kwargs = {}
+ style_kwargs.update(_user_figsize(args.width, args.height))
+ if args.fontsize:
+ style_kwargs.update({"font.size": args.fontsize})
+
+ style_list = get_style_list(
+ no_base_style=args.no_base_style, styles=args.style, **style_kwargs
+ )
+ matplotlib.style.use(style_list)
+
+ if args.sigma:
+ sigma = args.sigma
+ elif args.fwhm:
+ sigma = args.fwhm / (2 * sqrt(2 * log(2)))
+ else:
+ sigma = None
+
+ if args.action in ["automatic-plot"]:
+ describe.plot_cohps(
+ ylim=args.ylim,
+ xlim=args.xlim,
+ integrated=args.integrated,
+ save=args.save_plot is not None,
+ filename=args.save_plot,
+ title=args.title,
+ sigma=sigma,
+ hide=args.hideplot,
+ )
+
+ if args.action in [
+ "automatic-plot-ia",
+ "automaticplotia",
+ "auto-plot-ia",
+ "autoplotia",
+ ]:
+ describe.plot_interactive_cohps(
+ ylim=args.ylim,
+ xlim=args.xlim,
+ integrated=args.integrated,
+ save_as_html=args.save_plot,
+ filename=args.save_plot,
+ title=args.title,
+ sigma=sigma,
+ hide=args.hideplot,
+ label_resolved=args.labelresolved,
+ orbital_resolved=args.orbitalplot,
+ )
+
+ if args.action == "plot":
+ if args.cobis:
+ filename = args.cohpcar.parent / "COBICAR.lobster"
+ if not filename.exists():
+ filename = filename.with_name(filename.name + ".gz")
+ options = {"are_cobis": True, "are_coops": False}
+ elif args.coops:
+ filename = args.cohpcar.parent / "COOPCAR.lobster"
+ if not filename.exists():
+ filename = filename.with_name(filename.name + ".gz")
+ options = {"are_cobis": False, "are_coops": True}
+ else:
+ filename = args.cohpcar.parent / "COHPCAR.lobster"
+ if not filename.exists():
+ filename = filename.with_name(filename.name + ".gz")
+ options = {"are_cobis": False, "are_coops": False}
+
+ completecohp = CompleteCohp.from_file(
+ fmt="LOBSTER", filename=filename, structure_file=args.poscar, **options
+ )
+ cp = PlainCohpPlotter(**options)
+
+ if not args.summed:
+ # TODO: add checks for label in allowed labels -> print all labels
+ # TODO: add check if args.oribtalwise is exactly as long as labels
+ # TODO: add check if orbital is in args.orbitalwise
+
+ for label in args.bond_numbers:
+ if str(label) not in completecohp.bonds:
+ raise IndexError(
+ "The provided bond label "
+ + str(label)
+ + " is not available in ICO**LIST.lobster.\n "
+ "Allowed options are in this list: \n"
+ + str([int(listi) for listi in list(completecohp.bonds.keys())])
+ )
+
+ if not args.orbitalwise:
+ for label in args.bond_numbers:
+ cp.add_cohp(label, completecohp.get_cohp_by_label(label=str(label)))
+ else:
+ if len(args.bond_numbers) != len(args.orbitalwise):
+ raise IndexError(
+ "Please provide as mainy orbitals as bond labels,"
+ " e.g., lobsterpy plot 1 1 --orbitalwise '2s-2s' '2s-2px'"
+ )
+
+ for ilabel, label in enumerate(args.bond_numbers):
+ orbitals = args.orbitalwise[ilabel]
+
+ availableorbitals = list(
+ completecohp.orb_res_cohp[str(label)].keys()
+ )
+ orbitaloptions = [*availableorbitals, "all"]
+
+ if orbitals not in orbitaloptions:
+ raise IndexError(
+ "Orbital in not available for current bond. \n"
+ "For bond "
+ + str(label)
+ + " only the following orbital options are available: \n"
+ + str(orbitaloptions)
+ )
+
+ if orbitals != "all":
+ cp.add_cohp(
+ str(label) + ": " + orbitals,
+ completecohp.get_orbital_resolved_cohp(
+ label=str(label), orbitals=orbitals
+ ),
+ )
+ else:
+ for orbitals in availableorbitals:
+ cp.add_cohp(
+ str(label) + ": " + orbitals,
+ completecohp.get_orbital_resolved_cohp(
+ label=str(label), orbitals=orbitals
+ ),
+ )
+ else:
+ cp.add_cohp(
+ str(args.bond_numbers),
+ completecohp.get_summed_cohp_by_label_list(
+ label_list=[str(label) for label in args.bond_numbers]
+ ),
+ )
+
+ plt = cp.get_plot(
+ integrated=args.integrated, xlim=args.xlim, ylim=args.ylim, sigma=sigma
+ )
+
+ ax = plt.gca()
+ ax.set_title(args.title)
+
+ if not args.hideplot and not args.save_plot:
+ plt.show()
+ elif args.save_plot and not args.hideplot:
+ plt.show()
+ fig = plt.gcf()
+ fig.savefig(args.save_plot)
+ if args.save_plot and args.hideplot:
+ plt.savefig(args.save_plot)
+
+ if args.action in ["create-inputs", "createinputs"]:
+ from pymatgen.core.structure import Structure
+ from pymatgen.io.lobster import Lobsterin
+
+ # Check for .gz files exist for default values and update accordingly
+ default_files = {
+ "poscar": "POSCAR",
+ "potcar": "POTCAR",
+ "incar": "INCAR",
+ }
+
+ for arg_name in default_files:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "Files necessary for creating puts for LOBSTER calcs not found in the current directory."
+ )
+
+ if args.userbasis is None:
+ # This will rely on standard basis files as stored in pymatgen
+
+ potcar_names = Lobsterin._get_potcar_symbols(POTCAR_input=args.potcar)
+
+ list_basis_dict = Lobsterin.get_all_possible_basis_functions(
+ structure=Structure.from_file(args.poscar), potcar_symbols=potcar_names
+ )
+
+ for ibasis, basis_dict in enumerate(list_basis_dict):
+ lobsterinput = Lobsterin.standard_calculations_from_vasp_files(
+ args.poscar,
+ args.incar,
+ None,
+ option="standard",
+ dict_for_basis=basis_dict,
+ )
+
+ lobsterin_path = Path(str(args.lobsterinout) + "-" + str(ibasis))
+ incar_path = Path(str(args.incarout) + "-" + str(ibasis))
+
+ if (not lobsterin_path.is_file() and not incar_path.is_file()) or (
+ args.overwrite
+ ):
+ lobsterinput.write_lobsterin(lobsterin_path)
+ lobsterinput.write_INCAR(
+ incar_input=args.incar,
+ incar_output=incar_path,
+ poscar_input=args.poscar,
+ isym=0,
+ )
+ else:
+ raise ValueError(
+ 'please use "--overwrite" if you would like to overwrite existing lobster inputs'
+ )
+ else:
+ # convert list userbasis to dict
+ userbasis = {}
+ for userbasis_single in args.userbasis:
+ userbasis[userbasis_single[0]] = userbasis_single[1]
+
+ lobsterinput = Lobsterin.standard_calculations_from_vasp_files(
+ args.poscar,
+ args.incar,
+ None,
+ option="standard",
+ dict_for_basis=userbasis,
+ )
+
+ lobsterin_path = Path(str(args.lobsterinout) + "-" + str(0))
+ incar_path = Path(str(args.incarout) + "-" + str(0))
+
+ if (not lobsterin_path.is_file() and not incar_path.is_file()) or (
+ args.overwrite
+ ):
+ lobsterinput.write_lobsterin(lobsterin_path)
+ lobsterinput.write_INCAR(
+ incar_input=args.incar,
+ incar_output=incar_path,
+ poscar_input=args.poscar,
+ isym=0,
+ )
+ else:
+ raise ValueError(
+ 'please use "--overwrite" if you would like to overwrite existing lobster inputs'
+ )
+
+ if args.action in ["calc-description"]:
+ # Check for .gz files exist for default values and update accordingly
+ mandatory_files = {
+ "poscar": "POSCAR",
+ "lobsterin": "lobsterin",
+ "lobsterout": "lobsterout",
+ }
+
+ for arg_name in mandatory_files:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "Mandatory files necessary for LOBSTER calc quality not found in the current directory."
+ )
+
+ optional_file = {
+ "bandoverlaps": "bandOverlaps.lobster",
+ "potcar": "POTCAR",
+ }
+
+ for arg_name in optional_file:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+
+ bva_comp = args.bvacomp
+
+ if bva_comp:
+ bva_files = {
+ "charge": "CHARGE.lobster",
+ }
+ for arg_name in bva_files:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "BVA charge requested but CHARGE.lobster file not found."
+ )
+
+ dos_comparison = args.doscomp
+
+ if dos_comparison:
+ dos_files = {
+ "doscar": "DOSCAR.lobster",
+ "vasprun": "vasprun.xml",
+ }
+
+ for arg_name in dos_files:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ "DOS comparisons requested but DOSCAR.lobster, vasprun.xml file not found."
+ )
+ potcar_file_path = args.potcar
+
+ quality_dict = Analysis.get_lobster_calc_quality_summary(
+ path_to_poscar=args.poscar,
+ path_to_charge=args.charge,
+ path_to_lobsterout=args.lobsterout,
+ path_to_lobsterin=args.lobsterin,
+ path_to_potcar=None if not potcar_file_path.exists() else potcar_file_path,
+ potcar_symbols=args.potcarsymbols,
+ path_to_bandoverlaps=args.bandoverlaps,
+ dos_comparison=dos_comparison,
+ bva_comp=bva_comp,
+ path_to_doscar=args.doscar,
+ e_range=args.erange,
+ n_bins=args.nbins,
+ path_to_vasprun=args.vasprun,
+ )
+
+ quality_text = Description.get_calc_quality_description(quality_dict)
+ Description.write_calc_quality_description(quality_text)
+
+ if args.calcqualityjson is not None:
+ with open(args.calcqualityjson, "w") as fd:
+ json.dump(quality_dict, fd)
+
+ if args.action in ["plot-dos", "plotdos"]:
+ mandatory_files = {
+ "doscar": "DOSCAR.lobster",
+ "poscar": "POSCAR",
+ }
+
+ for arg_name in mandatory_files:
+ file_path = getattr(args, arg_name)
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ setattr(args, arg_name, gz_file_path)
+ else:
+ raise ValueError(
+ f"{file_path.name} necessary for plotting DOS not found in "
+ f"the current directory."
+ )
+
+ from pymatgen.io.lobster import Doscar
+
+ lobs_dos = Doscar(doscar=args.doscar, structure_file=args.poscar).completedos
+
+ dos_plotter = PlainDosPlotter(summed=args.summedspins, sigma=args.sigma)
+ if args.addtotaldos:
+ dos_plotter.add_dos(dos=lobs_dos, label="Total DOS")
+ if args.spddos:
+ dos_plotter.add_dos_dict(dos_dict=lobs_dos.get_spd_dos())
+
+ if args.elementdos:
+ dos_plotter.add_dos_dict(dos_dict=lobs_dos.get_element_dos())
+
+ if args.element:
+ for element in args.element:
+ element_spddos = lobs_dos.get_element_spd_dos(el=element)
+ for orbital, dos in element_spddos.items():
+ label = f"{element}: {orbital.name}"
+ dos_plotter.add_dos_dict(dos_dict={label: dos})
+
+ if args.site is not None and args.orbital:
+ if len(args.site) > len(args.orbital):
+ for site in args.site:
+ for orbital in args.orbital:
+ dos_plotter.add_site_orbital_dos(
+ site_index=site, orbital=orbital, dos=lobs_dos
+ )
+ elif len(args.orbital) > len(args.site):
+ for orbital in args.orbital:
+ for site in args.site:
+ dos_plotter.add_site_orbital_dos(
+ site_index=site, orbital=orbital, dos=lobs_dos
+ )
+ else:
+ for site, orbital in zip(args.site, args.orbital):
+ dos_plotter.add_site_orbital_dos(
+ site_index=site, orbital=orbital, dos=lobs_dos
+ )
+
+ elif (args.site is None or not args.orbital) and (
+ not args.element and not args.spddos and not args.elementdos
+ ):
+ raise ValueError(
+ "Please set both args i.e site and orbital to generate the plot"
+ )
+
+ plt = dos_plotter.get_plot(
+ xlim=args.xlim,
+ ylim=args.ylim,
+ beta_dashed=True,
+ invert_axes=args.invertaxis,
+ )
+
+ ax = plt.gca()
+ ax.set_title(args.title)
+
+ if not args.hideplot and not args.save_plot:
+ plt.show()
+ elif args.save_plot and not args.hideplot:
+ plt.show()
+ fig = plt.gcf()
+ fig.savefig(args.save_plot)
+ if args.save_plot and args.hideplot:
+ plt.savefig(args.save_plot)
+
+ if args.action in ["plot-icohps-distances", "ploticohpsdistances"]:
+ if args.cobis:
+ filename = args.icohplist.parent / "ICOBILIST.lobster"
+ if not filename.exists():
+ filename = filename.with_name(filename.name + ".gz")
+ options = {"are_cobis": True, "are_coops": False}
+ elif args.coops:
+ filename = args.icohplist.parent / "ICOOPLIST.lobster"
+ if not filename.exists():
+ filename = filename.with_name(filename.name + ".gz")
+ options = {"are_cobis": False, "are_coops": True}
+ else:
+ filename = args.icohplist.parent / "ICOHPLIST.lobster"
+ if not filename.exists():
+ filename = filename.with_name(filename.name + ".gz")
+ options = {"are_cobis": False, "are_coops": False}
+
+ icohpcollection = Icohplist(filename=filename, **options).icohpcollection
+ icohp_plotter = IcohpDistancePlotter(**options)
+
+ icohp_plotter.add_icohps(icohpcollection=icohpcollection, label="")
+
+ plt = icohp_plotter.get_plot(xlim=args.xlim, ylim=args.ylim)
+
+ ax = plt.gca()
+ ax.set_title(args.title)
+
+ if not args.hideplot and not args.save_plot:
+ plt.show()
+ elif args.save_plot and not args.hideplot:
+ plt.show()
+ fig = plt.gcf()
+ fig.savefig(args.save_plot)
+ if args.save_plot and args.hideplot:
+ plt.savefig(args.save_plot)
+
+
+if __name__ == "__main__":
+ main()
+
+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""This module defines classes to analyze the COHPs automatically."""
+from __future__ import annotations
+
+import warnings
+from collections import Counter
+from pathlib import Path
+
+import numpy as np
+from pymatgen.analysis.bond_valence import BVAnalyzer
+from pymatgen.core.structure import Structure
+from pymatgen.electronic_structure.cohp import CompleteCohp
+from pymatgen.electronic_structure.core import Spin
+from pymatgen.io.lobster import Bandoverlaps, Charge, Doscar, Lobsterin, Lobsterout
+from pymatgen.io.lobster.lobsterenv import LobsterNeighbors
+from pymatgen.io.vasp.outputs import Vasprun
+from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
+
+
+[docs]class Analysis:
+ """
+ Class to analyze COHP/COOP/COBI data from Lobster.
+
+ Attributes:
+ condensed_bonding_analysis: dict including a summary of the most important bonding properties
+ final_dict_bonds: dict including information on ICOHPs per bond type
+ final_dict_ions: dict including information on environments of cations
+ chemenv: pymatgen.io.lobster.lobsterenv.LobsterNeighbors object
+ lse: LightStructureEnvironment from pymatgen
+ cutoff_icohp: Cutoff in percentage for evaluating neighbors based on ICOHP values.
+ cutoff_icohp*max_icohp limits the number of considered environments
+ anion_types: Set of Element objects from pymatgen
+ list_equivalent_sites: list of site indices of sites that indicate which sites are equivalent
+ e.g., [0 1 2 2 2] where site 0, 1, 2 indicate sites that are independent from each other
+ path_to_charge: str that describes the path to CHARGE.lobster
+ path_to_cohpcar: str that describes the path to COHPCAR.lobster
+ path_to_icohplist: str that describes the path to ICOHPLIST.lobster
+ path_to_poscar: str that describes path to POSCAR
+ path_to_madelung: str that describes path to POSCAR
+ are_cobis : bool indicating if file contains COBI/ICOBI data
+ are_coops : bool indicating if file contains COOP/ICOOP data
+ noise_cutoff : float that sets the lower limit of icohps or icoops or icobis considered
+ seq_cohps: list of cohps
+ seq_coord_ions: list of co-ordination environment strings for each cation
+ seq_equivalent_sites: seq of inequivalent sites
+ seq_ineq_ions: seq of inequivalent cations/sites in the structure
+ seq_infos_bonds (list): information on cation anion bonds (lists
+ of pymatgen.io.lobster.lobsterenv.ICOHPNeighborsInfo)
+ spg: space group information
+ structure: Structure object
+ type_charge: which charges are considered here
+ orbital_resolved: bool indicating whether analysis is performed
+ orbital wise
+ which_bonds: which bonds will be considered in analysis
+ """
+
+ def __init__(
+ self,
+ path_to_poscar: str | Path,
+ path_to_icohplist: str | Path,
+ path_to_cohpcar: str | Path,
+ path_to_charge: str | Path | None = None,
+ path_to_madelung: str | Path | None = None,
+ which_bonds: str = "cation-anion",
+ cutoff_icohp: float = 0.1,
+ noise_cutoff: float = 0.1,
+ orbital_cutoff: float = 0.05,
+ summed_spins=True,
+ are_cobis=False,
+ are_coops=False,
+ type_charge=None,
+ start=None,
+ orbital_resolved=False,
+ ):
+ """
+ Automatically analyze bonding information with this class.
+
+ Args:
+ path_to_poscar: path to POSCAR (e.g., "POSCAR")
+ path_to_icohplist: path to ICOHPLIST.lobster (e.g., "ICOHPLIST.lobster")
+ path_to_cohpcar: path to COHPCAR.lobster (e.g., "COHPCAR.lobster")
+ path_to_charge: path to CHARGE.lobster (e.g., "CHARGE.lobster")
+ path_to_madelung: path to MadelungEnergies.lobster (e.g., "MadelungEnergies.lobster")
+ are_cobis : bool indicating if file contains COBI/ICOBI data
+ are_coops : bool indicating if file contains COOP/ICOOP data
+ noise_cutoff : float that sets the lower limit of icohps or icoops or icobis considered
+ orbital_cutoff : float that sets the minimum percentage for the orbital resolved analysis.
+ (Affects only when orbital_resolved argument is set to True) Set it to 0 to get results
+ of all orbitals in the detected relevant bonds. Default is to 0.05 i.e. only analyzes
+ if orbital contribution is 5 % or more.
+ which_bonds: selects which kind of bonds are analyzed. "cation-anion" is the default
+ cutoff_icohp: only bonds that are stronger than cutoff_icohp*strongest ICOHP will be considered
+ summed_spins: if true, spins will be summed
+ type_charge: If no path_to_charge is given, Valences will be used. Otherwise, Mulliken charges.
+ Löwdin charges cannot be selected at the moment.
+ orbital_resolved: bool indicating whether analysis is performed orbital wise
+ start: start energy for integration
+
+ """
+ self.start = start
+ self.path_to_poscar = path_to_poscar
+ self.path_to_icohplist = path_to_icohplist
+ self.path_to_cohpcar = path_to_cohpcar
+ self.which_bonds = which_bonds
+ self.cutoff_icohp = cutoff_icohp
+ self.orbital_cutoff = orbital_cutoff
+ self.path_to_charge = path_to_charge
+ self.path_to_madelung = path_to_madelung
+ self.are_cobis = are_cobis
+ self.are_coops = are_coops
+ self.noise_cutoff = noise_cutoff
+ self.setup_env()
+ self.get_information_all_bonds(summed_spins=summed_spins)
+ self.orbital_resolved = orbital_resolved
+
+ # This determines how cations and anions
+ if path_to_charge is None:
+ self.type_charge = "Valences"
+ else:
+ if type_charge is None:
+ self.type_charge = "Mulliken"
+ elif type_charge == "Mulliken":
+ self.type_charge = "Mulliken"
+ elif type_charge == "Löwdin":
+ raise ValueError(
+ "Only Mulliken charges can be used here at the moment. Implementation will follow."
+ )
+ else:
+ self.type_charge = "Valences"
+ print(
+ "type_charge cannot be read! Please use Mulliken/Löwdin. Now, we will use valences"
+ )
+
+ self.set_condensed_bonding_analysis()
+ self.set_summary_dicts()
+ self.path_to_madelung = path_to_madelung
+
+[docs] def setup_env(self):
+ """
+ Set up the light structure environments based on COHPs using this method.
+
+ Returns:
+ None
+
+ """
+ self.structure = Structure.from_file(self.path_to_poscar)
+ sga = SpacegroupAnalyzer(structure=self.structure)
+ symmetry_dataset = sga.get_symmetry_dataset()
+ equivalent_sites = symmetry_dataset["equivalent_atoms"]
+ self.list_equivalent_sites = equivalent_sites
+ self.seq_equivalent_sites = list(set(equivalent_sites))
+ self.spg = symmetry_dataset["international"]
+
+ if self.which_bonds == "cation-anion":
+ try:
+ self.chemenv = LobsterNeighbors(
+ filename_ICOHP=self.path_to_icohplist,
+ structure=Structure.from_file(self.path_to_poscar),
+ additional_condition=1,
+ perc_strength_ICOHP=self.cutoff_icohp,
+ filename_CHARGE=self.path_to_charge,
+ valences_from_charges=True,
+ adapt_extremum_to_add_cond=True,
+ are_cobis=self.are_cobis,
+ are_coops=self.are_coops,
+ noise_cutoff=self.noise_cutoff,
+ )
+ except ValueError as err:
+ if (
+ str(err) == "min() arg is an empty sequence"
+ or str(err)
+ == "All valences are equal to 0, additional_conditions 1, 3, 5 and 6 will not work"
+ ):
+ raise ValueError(
+ "Consider switching to an analysis of all bonds and not only cation-anion bonds."
+ " It looks like no cations are detected."
+ )
+ raise err
+ elif self.which_bonds == "all":
+ # raise ValueError("only cation anion bonds implemented so far")
+ self.chemenv = LobsterNeighbors(
+ filename_ICOHP=self.path_to_icohplist,
+ structure=Structure.from_file(self.path_to_poscar),
+ additional_condition=0,
+ perc_strength_ICOHP=self.cutoff_icohp,
+ filename_CHARGE=self.path_to_charge,
+ valences_from_charges=True,
+ adapt_extremum_to_add_cond=True,
+ are_cobis=self.are_cobis,
+ are_coops=self.are_coops,
+ noise_cutoff=self.noise_cutoff,
+ )
+
+ else:
+ raise ValueError("only cation anion and all bonds implemented so far")
+
+ # determine cations and anions
+ try:
+ if self.which_bonds == "cation-anion":
+ self.lse = self.chemenv.get_light_structure_environment(
+ only_cation_environments=True
+ )
+ elif self.which_bonds == "all":
+ self.lse = self.chemenv.get_light_structure_environment(
+ only_cation_environments=False
+ )
+ except ValueError:
+
+ class Lse:
+ """Test class when error was raised."""
+
+ def __init__(self, chemenv, valences=None):
+ """
+ Test class when error was raised.
+
+ Args:
+ chemenv (LobsterNeighbors): LobsterNeighbors object
+ valences: list of valences
+
+ """
+ if valences is None:
+ self.coordination_environments = [
+ [{"ce_symbol": str(len(coord))}] for coord in chemenv
+ ]
+ else:
+ self.coordination_environments = []
+
+ for val, coord in zip(valences, chemenv):
+ if val >= 0.0:
+ self.coordination_environments.append(
+ [{"ce_symbol": str(len(coord))}]
+ )
+ else:
+ self.coordination_environments.append(
+ [{"ce_symbol": None}]
+ )
+
+ if self.which_bonds == "all":
+ self.lse = Lse(self.chemenv.list_coords)
+ elif self.which_bonds == "cation-anion":
+ # make a new list
+ self.lse = Lse(self.chemenv.list_coords, self.chemenv.valences)
+
+[docs] def get_information_all_bonds(self, summed_spins=True):
+ """
+ Gather all information on the bonds within the compound with this method.
+
+ Returns:
+ None
+
+ """
+ if self.which_bonds == "cation-anion":
+ # this will only analyze cation anion bonds which simplifies the analysis
+ self.seq_ineq_ions = []
+ self.seq_coord_ions = []
+ self.seq_infos_bonds = []
+ self.seq_labels_cohps = []
+ self.seq_cohps = []
+ # only_bonds_to
+
+ self.anion_types = self.chemenv.get_anion_types()
+ for ice, ce in enumerate(self.lse.coordination_environments):
+ # only look at inequivalent sites (use of symmetry to speed everything up!)!
+ # only look at those cations that have cation-anion bonds
+ if ice in self.seq_equivalent_sites and ce[0]["ce_symbol"] is not None:
+ self.seq_ineq_ions.append(ice)
+
+ self.seq_coord_ions.append(ce[0]["ce_symbol"])
+ self.seq_infos_bonds.append(
+ self.chemenv.get_info_icohps_to_neighbors([ice])
+ )
+
+ aniontype_labels = []
+ aniontype_cohps = []
+
+ # go through all anions in the structure!
+ for anion in self.anion_types:
+ # get labels and summed cohp objects
+ labels, summedcohps = self.chemenv.get_info_cohps_to_neighbors(
+ self.path_to_cohpcar,
+ [ice],
+ summed_spin_channels=summed_spins,
+ per_bond=False,
+ only_bonds_to=[str(anion)],
+ )
+
+ aniontype_labels.append(labels)
+ aniontype_cohps.append(summedcohps)
+
+ self.seq_labels_cohps.append(aniontype_labels)
+ self.seq_cohps.append(aniontype_cohps)
+
+ elif self.which_bonds == "all":
+ # this will only analyze all bonds
+
+ self.seq_ineq_ions = []
+ self.seq_coord_ions = []
+ self.seq_infos_bonds = []
+ self.seq_labels_cohps = []
+ self.seq_cohps = []
+ # only_bonds_to
+ self.elements = self.structure.composition.elements
+ # self.anion_types = self.chemenv.get_anion_types()
+ for ice, ce in enumerate(self.lse.coordination_environments):
+ # only look at inequivalent sites (use of symmetry to speed everything up!)!
+ # only look at those cations that have cation-anion bonds
+ if ice in self.seq_equivalent_sites and ce[0]["ce_symbol"] is not None:
+ self.seq_ineq_ions.append(ice)
+ self.seq_coord_ions.append(ce[0]["ce_symbol"])
+ self.seq_infos_bonds.append(
+ self.chemenv.get_info_icohps_to_neighbors([ice])
+ )
+
+ type_labels = []
+ type_cohps = []
+
+ for element in self.elements:
+ # get labels and summed cohp objects
+ labels, summedcohps = self.chemenv.get_info_cohps_to_neighbors(
+ self.path_to_cohpcar,
+ [ice],
+ onlycation_isites=False,
+ summed_spin_channels=summed_spins,
+ per_bond=False,
+ only_bonds_to=[str(element)],
+ )
+
+ type_labels.append(labels)
+ type_cohps.append(summedcohps)
+
+ self.seq_labels_cohps.append(type_labels)
+ self.seq_cohps.append(type_cohps)
+
+[docs] def get_site_bond_resolved_labels(self):
+ """
+ Return relevant bond labels for each symmetrically independent site.
+
+ Returns:
+ dict with bond labels for each site, e.g.
+ {'Na1: Na-Cl': ['21', '23', '24', '27', '28', '30']}
+
+ """
+ bonds = [[] for _ in range(len(self.seq_infos_bonds))] # type: ignore
+ labels = [[] for _ in range(len(self.seq_infos_bonds))] # type: ignore
+ for inx, bond_info in enumerate(self.seq_infos_bonds):
+ for ixx, val in enumerate(bond_info.atoms):
+ label_srt = sorted(val.copy())
+ bonds[inx].append(
+ self.structure.sites[bond_info.central_isites[0]].species_string
+ + str(bond_info.central_isites[0] + 1)
+ + ": "
+ + label_srt[0].strip("0123456789")
+ + "-"
+ + label_srt[1].strip("0123456789")
+ )
+ labels[inx].append(bond_info.labels[ixx])
+
+ label_data = {}
+ for indx, atom_pairs in enumerate(bonds):
+ searched_atom_pairs = set(atom_pairs)
+ for search_item in searched_atom_pairs:
+ indices = [
+ i for i, pair in enumerate(atom_pairs) if pair == search_item
+ ]
+ filtered_bond_label_list = [labels[indx][i] for i in indices]
+ label_data.update({search_item: filtered_bond_label_list})
+
+ return label_data
+
+[docs] def _get_orbital_resolved_data(
+ self, nameion, iion, labels, bond_resolved_labels, type_pop
+ ):
+ """
+ Retrieve orbital-wise analysis data.
+
+ Args:
+ nameion: name of symmetrically relevant cation or anion
+ iion: index of symmetrically relevant cation or anion
+ labels: list of bond label names
+ bond_resolved_labels: dict of bond labels from ICOHPLIST resolved for each bond
+ type_pop: population type analyzed. e.g. COHP or COOP or COBI
+
+ Returns:
+ dict consisting of relevant orbitals (contribution > 5 % to overall ICOHP or ICOBI or ICOOP),
+ bonding and antibonding percentages with bond label names as keys.
+ """
+ orb_resolved_bond_info = {}
+ for label in labels:
+ if label is not None:
+ bond_resolved_label_key = (
+ nameion + str(iion + 1) + ":" + label.split("x")[-1]
+ )
+ bond_labels = bond_resolved_labels[bond_resolved_label_key]
+ available_orbitals = list(
+ self.chemenv.completecohp.orb_res_cohp[bond_labels[0]].keys()
+ )
+ # initialize empty list to store orb paris for bonding,
+ # antibonding integrals and percentages
+ bndg_orb_pair_list = []
+ bndg_orb_integral_list = []
+ bndg_orb_perc_list = []
+ bndg_orb_icohp_list = []
+ antibndg_orb_integral_list = []
+ antibndg_orb_perc_list = []
+ antibndg_orb_pair_list = []
+ antibndg_orb_icohp_list = []
+
+ # get total summed cohps using label list
+ cohp_summed = self.chemenv.completecohp.get_summed_cohp_by_label_list(
+ label_list=bond_labels
+ )
+ if type_pop.lower() == "cohp":
+ (
+ antibndg_tot,
+ per_anti_tot,
+ bndg_tot,
+ per_bndg_tot,
+ ) = self._integrate_antbdstates_below_efermi(
+ cohp=cohp_summed, start=self.start
+ )
+ else:
+ (
+ bndg_tot,
+ per_bndg_tot,
+ antibndg_tot,
+ per_anti_tot,
+ ) = self._integrate_antbdstates_below_efermi(
+ cohp=cohp_summed, start=self.start
+ )
+
+ orb_bonding_dict_data = {}
+ # For each orbital collect the contributions of summed bonding
+ # and antibonding interactions separately
+ for orb in available_orbitals:
+ cohp_summed_orb = self.chemenv.completecohp.get_summed_cohp_by_label_and_orbital_list(
+ label_list=bond_labels, orbital_list=[orb] * len(bond_labels)
+ )
+
+ if type_pop.lower() == "cohp":
+ (
+ antibndg_orb,
+ per_anti_orb,
+ bndg_orb,
+ per_bndg_orb,
+ ) = self._integrate_antbdstates_below_efermi(
+ cohp=cohp_summed_orb, start=self.start
+ )
+ else:
+ (
+ bndg_orb,
+ per_bndg_orb,
+ antibndg_orb,
+ per_anti_orb,
+ ) = self._integrate_antbdstates_below_efermi(
+ cohp=cohp_summed_orb, start=self.start
+ )
+
+ # replace nan values with zero (tackle numerical integration issues)
+ bndg_orb = bndg_orb if not np.isnan(bndg_orb) else 0
+ per_bndg_orb = per_bndg_orb if not np.isnan(per_bndg_orb) else 0
+ bndg_tot = bndg_tot if not np.isnan(bndg_tot) else 0
+ per_bndg_tot = per_bndg_tot if not np.isnan(per_bndg_tot) else 0
+ # skip collecting orb contributions if no summed bonding contribution exists
+ if bndg_tot > 0:
+ orb_icohps_bndg = []
+ for bond_label in bond_labels:
+ orb_icohp_bn = (
+ self.chemenv.Icohpcollection.get_icohp_by_label(
+ label=bond_label, orbitals=orb
+ )
+ )
+ orb_icohps_bndg.append(orb_icohp_bn)
+ bndg_orb_pair_list.append(orb)
+ bndg_orb_icohp_list.append(orb_icohps_bndg)
+ bndg_orb_integral_list.append(bndg_orb)
+ bndg_orb_perc_list.append(per_bndg_orb)
+
+ # replace nan values with zero (tackle numerical integration issues)
+ antibndg_orb = antibndg_orb if not np.isnan(antibndg_orb) else 0
+ per_anti_orb = per_anti_orb if not np.isnan(per_anti_orb) else 0
+ antibndg_tot = antibndg_tot if not np.isnan(antibndg_tot) else 0
+ per_anti_tot = per_anti_tot if not np.isnan(per_anti_tot) else 0
+ # skip collecting orb contributions if no summed antibonding contribution exists
+ if antibndg_tot > 0:
+ orb_icohps_anti = []
+ for bond_label in bond_labels:
+ orb_icohp_an = (
+ self.chemenv.Icohpcollection.get_icohp_by_label(
+ label=bond_label, orbitals=orb
+ )
+ )
+ orb_icohps_anti.append(orb_icohp_an)
+
+ antibndg_orb_pair_list.append(orb)
+ antibndg_orb_icohp_list.append(orb_icohps_anti)
+ antibndg_orb_integral_list.append(antibndg_orb)
+ antibndg_orb_perc_list.append(per_anti_orb)
+
+ # Populate the dictionary with relevant orbitals for bonding interactions
+ for inx, bndg_orb_pair in enumerate(bndg_orb_pair_list):
+ bndg_contri_perc = round(
+ bndg_orb_integral_list[inx] / sum(bndg_orb_integral_list), 2
+ )
+ # filter out very small bonding interactions (<self.orbital_cutoff)
+ if bndg_contri_perc > self.orbital_cutoff:
+ if bndg_orb_pair in orb_bonding_dict_data:
+ orb_bonding_dict_data[bndg_orb_pair].update(
+ {
+ "orb_contribution_perc_bonding": bndg_contri_perc,
+ "bonding": {
+ "integral": bndg_orb_integral_list[inx],
+ "perc": bndg_orb_perc_list[inx],
+ },
+ }
+ )
+ else:
+ orb_bonding_dict_data[bndg_orb_pair] = {
+ f"I{type_pop}_mean": round(
+ np.mean(bndg_orb_icohp_list[inx]), 4
+ ),
+ f"I{type_pop}_sum": round(
+ np.sum(bndg_orb_icohp_list[inx]), 4
+ ),
+ "orb_contribution_perc_bonding": round(
+ bndg_orb_integral_list[inx]
+ / sum(bndg_orb_integral_list),
+ 2,
+ ),
+ "bonding": {
+ "integral": bndg_orb_integral_list[inx],
+ "perc": bndg_orb_perc_list[inx],
+ },
+ }
+
+ # Populate the dictionary with relevant orbitals for antibonding interactions
+ for inx, antibndg_orb_pair in enumerate(antibndg_orb_pair_list):
+ antibndg_contri_perc = round(
+ antibndg_orb_integral_list[inx]
+ / sum(antibndg_orb_integral_list),
+ 2,
+ )
+ # filter out very small antibonding interactions (<self.orbital_cutoff)
+ if antibndg_contri_perc > self.orbital_cutoff:
+ if antibndg_orb_pair in orb_bonding_dict_data:
+ orb_bonding_dict_data[antibndg_orb_pair].update(
+ {
+ "orb_contribution_perc_antibonding": round(
+ antibndg_orb_integral_list[inx]
+ / sum(antibndg_orb_integral_list),
+ 2,
+ ),
+ "antibonding": {
+ "integral": antibndg_orb_integral_list[inx],
+ "perc": antibndg_orb_perc_list[inx],
+ },
+ }
+ )
+ else:
+ orb_bonding_dict_data[antibndg_orb_pair] = {
+ f"I{type_pop}_mean": round(
+ np.mean(antibndg_orb_icohp_list[inx]), 4
+ ),
+ f"I{type_pop}_sum": round(
+ np.sum(antibndg_orb_icohp_list[inx]), 4
+ ),
+ "orb_contribution_perc_antibonding": round(
+ antibndg_orb_integral_list[inx]
+ / sum(antibndg_orb_integral_list),
+ 2,
+ ),
+ "antibonding": {
+ "integral": antibndg_orb_integral_list[inx],
+ "perc": antibndg_orb_perc_list[inx],
+ },
+ }
+
+ orb_bonding_dict_data["relevant_bonds"] = bond_labels
+
+ orb_resolved_bond_info[bond_resolved_label_key] = orb_bonding_dict_data
+
+ return orb_resolved_bond_info
+
+[docs] def _get_bond_resolved_data_stats(self, orb_resolved_bond_data: dict):
+ """
+ Retrieve the maximum bonding and anti-bonding orbital contributions.
+
+ Args:
+ orb_resolved_bond_data: A dictionary with orbital names as keys and corresponding bonding data
+
+ Returns:
+ dict with orbital data stats the site for relevant orbitals, e.g.
+ {'orbital_summary_stats': {'max_bonding_contribution': {'2s-3s': 0.68},
+ 'max_antibonding_contribution': {'2s-2pz': 0.36}}}
+
+ """
+ # get max orbital bonding and contribution for the site
+ orb_pairs_bndg = []
+ orb_pairs_antibndg = []
+ orb_contri_bndg = []
+ orb_contri_antibndg = []
+ orbital_summary_stats = {"orbital_summary_stats": {}} # type: ignore
+ if orb_resolved_bond_data:
+ for orb_pair, data in orb_resolved_bond_data.items():
+ if "orb_contribution_perc_bonding" in data:
+ orb_pairs_bndg.append(orb_pair)
+ orb_contri_bndg.append(data["orb_contribution_perc_bonding"])
+
+ if "orb_contribution_perc_antibonding" in data:
+ orb_pairs_antibndg.append(orb_pair)
+ orb_contri_antibndg.append(
+ data["orb_contribution_perc_antibonding"]
+ )
+
+ if orb_contri_bndg:
+ max_orb_contri_bndg = max(orb_contri_bndg)
+ max_orb_contri_bndg_inxs = [
+ inx
+ for inx, orb_contri in enumerate(orb_contri_bndg)
+ if orb_contri == max_orb_contri_bndg
+ ]
+ max_orb_contri_bndg_dict = {}
+ for inx in max_orb_contri_bndg_inxs:
+ max_orb_contri_bndg_dict[orb_pairs_bndg[inx]] = orb_contri_bndg[inx]
+ orbital_summary_stats["orbital_summary_stats"][
+ "max_bonding_contribution"
+ ] = max_orb_contri_bndg_dict
+ if orb_contri_antibndg:
+ max_orb_contri_antibndg = max(orb_contri_antibndg)
+ max_antibndg_contri_inxs = [
+ inx
+ for inx, orb_anti_per in enumerate(orb_contri_antibndg)
+ if orb_anti_per == max_orb_contri_antibndg
+ ]
+ max_antibndg_contri_dict = {}
+ for inx in max_antibndg_contri_inxs:
+ max_antibndg_contri_dict[
+ orb_pairs_antibndg[inx]
+ ] = orb_contri_antibndg[inx]
+ orbital_summary_stats["orbital_summary_stats"][
+ "max_antibonding_contribution"
+ ] = max_antibndg_contri_dict
+
+ return orbital_summary_stats
+
+[docs] def get_site_orbital_resolved_labels(self):
+ """
+ Return relevant orbitals and bond labels for each symmetrically independent site.
+
+ Returns:
+ dict with bond labels for each site for relevant orbitals, e.g.
+ {'Na1: Na-Cl': {'3s-3s': ['21', '23', '24', '27', '28', '30']}
+
+ """
+ site_bond_labels = self.get_site_bond_resolved_labels()
+ orb_plot_data = {atom_pair: {} for atom_pair in site_bond_labels}
+ if self.orbital_resolved:
+ for site_index, cba_data in self.condensed_bonding_analysis[
+ "sites"
+ ].items():
+ for atom in cba_data["bonds"]:
+ for orb_pair in cba_data["bonds"][atom]["orbital_data"]:
+ if orb_pair not in ("orbital_summary_stats", "relevant_bonds"):
+ atom_pair = [cba_data["ion"], atom]
+ atom_pair.sort()
+ key = (
+ self.structure.sites[site_index].species_string
+ + str(site_index + 1)
+ + ": "
+ + "-".join(atom_pair)
+ )
+ label_list = site_bond_labels[key]
+ orb_plot_data[key].update({orb_pair: label_list})
+ else:
+ print(
+ "Please set orbital_resolved to True when instantiating Analysis object, "
+ "to get this data"
+ )
+
+ return orb_plot_data
+
+[docs] @staticmethod
+ def _get_strenghts_for_each_bond(pairs, strengths, nameion=None):
+ """
+ Return a dictionary of bond strengths.
+
+ Args:
+ pairs: list of list including labels for the atoms, e.g., [['O3', 'Cu1'], ['O3', 'Cu1']]
+ strengths (list of float): list that gives the icohp strengths as a float, [-1.86287, -1.86288]
+ nameion: string including the name of the cation in the list, e.g Cu1
+
+ Returns:
+ dict including inormation on icohps for each bond type, e.g.
+ {'Yb-Sb': [-1.59769, -2.14723, -1.7925, -1.60773, -1.80149, -2.14335]}
+
+
+ """
+ dict_strenghts = {}
+
+ for pair, strength in zip(pairs, strengths):
+ if nameion is not None:
+ new = [
+ LobsterNeighbors._split_string(pair[0])[0],
+ LobsterNeighbors._split_string(pair[1])[0],
+ ]
+ new = Analysis._sort_name(new, nameion)
+ string_here = new[0] + "-" + new[1]
+ else:
+ new = sorted(
+ [
+ LobsterNeighbors._split_string(pair[0])[0],
+ LobsterNeighbors._split_string(pair[1])[0],
+ ]
+ )
+ string_here = new[0] + "-" + new[1]
+
+ if string_here not in dict_strenghts:
+ dict_strenghts[string_here] = []
+ dict_strenghts[string_here].append(strength)
+ return dict_strenghts
+
+[docs] @staticmethod
+ def _sort_name(pair, nameion=None):
+ """
+ Place the cation first in a list of name strings.
+
+ Args:
+ pair: ["O","Cu"]
+ nameion: "Cu"
+
+ Returns:
+ will return list of str, e.g. ["Cu", "O"]
+
+ """
+ if nameion is not None:
+ new = []
+ if pair[0] == nameion:
+ new.append(pair[0])
+ new.append(pair[1])
+
+ elif pair[1] == nameion:
+ new.append(pair[1])
+ new.append(pair[0])
+
+ return new
+
+[docs] @staticmethod
+ def _sort_orbital_atom_pair(
+ atom_pair: list,
+ label: str,
+ complete_cohp: CompleteCohp,
+ orb_pair: str,
+ ):
+ """
+ Place the cation first in a list of name strings and add the associated orbital name alongside the atom name.
+
+ Args:
+ atom_pair: list of atom pair with cation first eg., ["Cl","Na"]
+ label: LOBSTER relevant bond label eg ., "3"
+ complete_cohp: pymatgen CompleteCohp object
+ orb_pair: relevant orbital pair eg., "2px-3s"
+
+ Returns:
+ will return list of str, e.g. ["Na(2px)", "Cl(3s)"]
+
+ """
+ orb_atom = {} # type: ignore
+ orb_pair_list = orb_pair.split("-")
+ # get orbital associated to the atom and store in a dict
+ for _inx, (site, site_orb) in enumerate(
+ zip(complete_cohp.bonds[label]["sites"], orb_pair_list)
+ ):
+ if (
+ site.species_string in orb_atom
+ ): # check necessary for bonds between same atoms
+ orb_atom[site.species_string].append(site_orb)
+ else:
+ orb_atom[site.species_string] = [site_orb]
+
+ orb_atom_list = []
+ # add orbital name next to atom_pair
+ for inx, atom in enumerate(atom_pair):
+ # check to ensure getting 2nd orbital if bond is between same atomic species
+ if inx == 1 and len(orb_atom.get(atom)) > 1: # type: ignore
+ atom_with_orb_name = f"{atom}({orb_atom.get(atom)[1]})" # type: ignore
+ else:
+ atom_with_orb_name = f"{atom}({orb_atom.get(atom)[0]})" # type: ignore
+ orb_atom_list.append(atom_with_orb_name)
+
+ return orb_atom_list
+
+[docs] def _get_antibdg_states(self, cohps, labels, nameion=None, limit=0.01):
+ """
+ Return a dictionary containing information on anti-bonding states.
+
+ e.g., similar to: {'Cu-O': True, 'Cu-F': True}
+
+ Args:
+ cohps: list of pymatgen.electronic_structure.cohp.Cohp objects
+ labels: ['2 x Cu-O', '4 x Cu-F']
+ nameion: string of the cation name, e.g. "Cu"
+ limit: limit to detect antibonding states
+
+ Returns:
+ dict including in formation on whether antibonding interactions exist,
+ e.g., {'Cu-O': True, 'Cu-F': True}
+
+
+ """
+ dict_antibd = {}
+ for label, cohp in zip(labels, cohps):
+ if label is not None:
+ if nameion is not None:
+ new = label.split(" ")[2].split("-")
+ sorted_new = self._sort_name(new, nameion)
+ new_label = sorted_new[0] + "-" + sorted_new[1]
+ else:
+ new = label.split(" ")[2].split("-")
+ sorted_new = sorted(new.copy())
+ new_label = sorted_new[0] + "-" + sorted_new[1]
+
+ antbd = cohp.has_antibnd_states_below_efermi(limit=limit)
+ if Spin.down in antbd:
+ dict_antibd[new_label] = antbd[Spin.up] or antbd[Spin.down]
+ else:
+ dict_antibd[new_label] = antbd[Spin.up]
+
+ return dict_antibd
+
+[docs] def _integrate_antbdstates_below_efermi_for_set_cohps(self, labels, cohps, nameion):
+ """
+ Return a dictionary containing information on antibonding states.
+
+ .. warning:: NEEDS MORE TESTS
+
+ It is important to note that only the energy range that has been computed can be considered
+ (i.e., this might not be all)
+
+ e.g. output: {'Cu-O': {'integral': 4.24374775705, 'perc': 5.7437713186999995},
+ 'Cu-F': {'integral': 3.07098300965, 'perc': 4.25800841445}}
+
+ Args:
+ cohps: list of pymatgen.electronic_structure.cohp.Cohp objects
+ labels: ['2 x Cu-O', '4 x Cu-F']
+ nameion: string of the cation name, e.g. "Cu"
+
+ Returns:
+ dict including in formation on whether antibonding interactions exist,
+ e.g., {'Cu-O': {'integral': 4.24374775705, 'perc': 5.7437713186999995},
+ 'Cu-F': {'integral': 3.07098300965, 'perc': 4.25800841445}}}
+ """
+ dict_bd_antibd = {}
+ for label, cohp in zip(labels, cohps):
+ if label is not None:
+ new = label.split(" ")[2].split("-")
+ sorted_new = self._sort_name(new, nameion)
+ new_label = sorted_new[0] + "-" + sorted_new[1]
+ if not self.are_cobis and not self.are_coops:
+ (
+ integral,
+ perc,
+ integral2,
+ perc2,
+ ) = self._integrate_antbdstates_below_efermi(cohp, start=self.start)
+ else:
+ (
+ integral2,
+ perc2,
+ integral,
+ perc,
+ ) = self._integrate_antbdstates_below_efermi(cohp, start=self.start)
+
+ if integral == 0 and integral2 != 0.0:
+ dict_bd_antibd[new_label] = {
+ "bonding": {"integral": integral2, "perc": perc2},
+ "antibonding": {"integral": integral, "perc": 0.0},
+ }
+ elif integral2 == 0.0 and integral != 0.0:
+ dict_bd_antibd[new_label] = {
+ "bonding": {"integral": integral2, "perc": 0.0},
+ "antibonding": {"integral": integral, "perc": perc},
+ }
+ elif integral == 0.0 and integral2 == 0.0:
+ dict_bd_antibd[new_label] = {
+ "bonding": {"integral": integral2, "perc": 0.0},
+ "antibonding": {"integral": integral, "perc": 0.0},
+ }
+ else:
+ dict_bd_antibd[new_label] = {
+ "bonding": {"integral": integral2, "perc": perc2},
+ "antibonding": {"integral": integral, "perc": perc},
+ }
+
+ return dict_bd_antibd
+
+[docs] def _integrate_antbdstates_below_efermi(self, cohp, start):
+ """
+ Integrate the cohp data to compute bonding and anti-bonding contribution below efermi.
+
+ .. warning:: NEEDS MORE TESTS
+
+ This integrates the whole COHP curve that has been computed.
+ The energy range is very important.
+ At present the energy range considered is dependent on COHPstartEnergy
+ set during lobster runs. The bonding / antibonding integral values are sensitive to this parameter.
+ If COHPstartEnergy value does not cover entire range of VASP calculations then
+ absolute value of ICOHP_sum might not be equivalent to (bonding- antibonding) integral values.
+
+ Args:
+ cohp: cohp object
+ start: integration start energy in eV , eg start = -15
+
+ Returns:
+ absolute value of antibonding, percentage value of antibonding,
+ absolute value of bonding and percentage value of bonding interactions
+ """
+ warnings.warn(
+ "The bonding, antibonding integral/percent values are numerical estimate."
+ " These values are sensitive to COHPstartEnergy parameter."
+ " If COHPstartEnergy value does not cover entire range of VASP calculations then"
+ " absolute value of ICOHP_sum might not be equivalent to (bonding- antibonding) integral values."
+ )
+
+ from scipy.integrate import trapezoid
+
+ def integrate_positive(y, x):
+ """
+ Integrate only bonding interactions of COHPs.
+
+ Args:
+ y: COHP values
+ x: Energy values
+
+ Returns:
+ integrated value of bonding interactions
+ """
+ y = np.asanyarray(y)
+ x = np.asanyarray(x)
+
+ bonding = trapezoid(y, x)
+
+ return np.round(bonding, 2)
+
+ def integrate_negative(y, x):
+ """
+ Integrate only anti-bonding interactions of COHPs.
+
+ Args:
+ y: COHP values
+ x: Energy values
+
+ Returns:
+ integrated value of anti-bonding interactions
+ """
+ y = np.asanyarray(y)
+ x = np.asanyarray(x)
+ antibonding = trapezoid(y, x)
+
+ return np.round(antibonding, 2)
+
+ # will integrate spin.up and spin.down only below efermi
+ energies_corrected = cohp.energies - cohp.efermi
+ summedcohp = (
+ cohp.cohp[Spin.up] + cohp.cohp[Spin.down]
+ if Spin.down in cohp.cohp
+ else cohp.cohp[Spin.up]
+ )
+
+ cohp_bf = []
+ en_bf = []
+
+ for i, en in enumerate(energies_corrected):
+ if (start is None) and en <= 0:
+ en_bf.append(en)
+ cohp_bf.append(-1 * summedcohp[i])
+ if (start is not None) and 0 >= en >= start:
+ en_bf.append(en)
+ cohp_bf.append(-1 * summedcohp[i])
+
+ # Separate the bonding and antibonding COHP values in separate lists
+ pos = []
+ en_pos = []
+ neg = []
+ en_neg = []
+
+ for i, scohp in enumerate(cohp_bf):
+ if scohp >= 0:
+ pos.append(scohp)
+ en_pos.append(energies_corrected[i])
+
+ for i, scohp in enumerate(cohp_bf):
+ if scohp <= 0:
+ neg.append(-1 * scohp)
+ en_neg.append(energies_corrected[i])
+
+ antibonding = integrate_negative(y=neg, x=en_neg)
+
+ bonding = integrate_positive(y=pos, x=en_pos)
+
+ return (
+ antibonding,
+ np.round(abs(antibonding) / (abs(bonding) + abs(antibonding)), 5),
+ bonding,
+ np.round(abs(bonding) / (abs(bonding) + abs(antibonding)), 5),
+ )
+
+[docs] def _get_pop_type(self):
+ """
+ Return the type of the input population file.
+
+ Returns:
+ A String of analysed population can be COOP/COBI/COHP
+ """
+ if self.are_cobis:
+ type_pop = "COBI"
+ elif self.are_coops:
+ type_pop = "COOP"
+ else:
+ type_pop = "COHP"
+
+ return type_pop
+
+[docs] @staticmethod
+ def _get_bond_dict(
+ bond_strength_dict,
+ small_antbd_dict,
+ nameion=None,
+ large_antbd_dict=None,
+ type_pop=None,
+ ):
+ """
+ Return a bond_dict that contains information for each site.
+
+ Args:
+ bond_strength_dict (dict): dict with bond names as key and lists of bond strengths as items
+ small_antbd_dict (dict): dict including if there are antibonding interactions, {'Yb-Sb': False}
+ nameion (str): name of the cation, e.g. Yb
+ large_antbd_dict: will be implemented later
+ type_pop: population type analyzed. eg. COHP
+
+ Returns:
+ Eg., if type_pop == 'COHP', will return
+ dict including information on the anion (as label) and the ICOHPs in the item of the dict
+ ICOHP_mean refers to the mean ICOHP in eV
+ ICOHP_sum refers to the sum of the ICOHPs in eV
+ has_antibdg_states_below_Efermi is True if there are antibonding interactions below Efermi
+ "number_of_bonds" will count the numbers of bonds to the cation
+
+ Example:
+ {'Sb': {'ICOHP_mean': '-1.85', 'ICOHP_sum': '-11.09',
+ 'has_antibdg_states_below_Efermi': False, 'number_of_bonds': 6}}
+
+ """
+ bond_dict = {}
+
+ for key, item in bond_strength_dict.items():
+ if nameion is not None:
+ a = key.split("-")[0]
+ b = key.split("-")[1]
+ if a == nameion:
+ key_here = b
+ elif b == nameion:
+ key_here = a
+
+ if large_antbd_dict is None:
+ bond_dict[key_here] = {
+ f"I{type_pop}_mean": str(round(np.mean(item), 2)),
+ f"I{type_pop}_sum": str(round(np.sum(item), 2)),
+ "has_antibdg_states_below_Efermi": small_antbd_dict[key],
+ "number_of_bonds": len(item),
+ }
+ else:
+ bond_dict[key_here] = {
+ f"I{type_pop}_mean": str(round(np.mean(item), 2)),
+ f"I{type_pop}_sum": str(round(np.sum(item), 2)),
+ "has_antibdg_states_below_Efermi": small_antbd_dict[key],
+ "number_of_bonds": len(item),
+ "perc_antibdg_states_below_Efermi": large_antbd_dict[key],
+ }
+
+ return bond_dict
+
+[docs] def set_condensed_bonding_analysis(self):
+ """
+ Condense the bonding analysis into a summary dictionary.
+
+ Returns:
+ None
+
+ """
+ self.condensed_bonding_analysis = {}
+ # which icohps are considered
+ if self.which_bonds == "cation-anion":
+ limit_icohps = self.chemenv._get_limit_from_extremum(
+ self.chemenv.Icohpcollection,
+ self.cutoff_icohp,
+ adapt_extremum_to_add_cond=True,
+ additional_condition=1,
+ )
+ elif self.which_bonds == "all":
+ limit_icohps = self.chemenv._get_limit_from_extremum(
+ self.chemenv.Icohpcollection,
+ self.cutoff_icohp,
+ adapt_extremum_to_add_cond=True,
+ additional_condition=0,
+ )
+ # formula of the compound
+ formula = str(self.structure.composition.reduced_formula)
+ # set population type
+ type_pop = self._get_pop_type()
+ # how many inequivalent cations are in the structure
+ if self.which_bonds == "cation-anion":
+ number_considered_ions = len(self.seq_ineq_ions)
+ elif self.which_bonds == "all":
+ number_considered_ions = len(self.seq_ineq_ions)
+
+ # what was the maximum bond lengths that was considered
+ max_bond_lengths = max(self.chemenv.Icohpcollection._list_length)
+
+ # what are the charges for the cations in the structure
+ charge_list = self.chemenv.valences
+
+ # dictionary including bonding information for each site
+ site_dict = {}
+ if self.which_bonds == "cation-anion":
+ for ication, ce, cation_anion_infos, labels, cohps in zip(
+ self.seq_ineq_ions,
+ self.seq_coord_ions,
+ self.seq_infos_bonds,
+ self.seq_labels_cohps,
+ self.seq_cohps,
+ ):
+ namecation = str(self.structure[ication].specie)
+
+ # This will compute the mean strengths of ICOHPs
+ mean_icohps = self._get_strenghts_for_each_bond(
+ pairs=cation_anion_infos[4],
+ strengths=cation_anion_infos[1],
+ nameion=namecation,
+ )
+ # pairs, strengths, nameion
+ # will collect if there are antibonding states present
+ antbdg = self._get_antibdg_states(cohps, labels, namecation)
+ dict_antibonding = (
+ self._integrate_antbdstates_below_efermi_for_set_cohps(
+ labels, cohps, nameion=namecation
+ )
+ )
+ bond_dict = self._get_bond_dict(
+ mean_icohps, antbdg, namecation, type_pop=type_pop
+ )
+ bond_resolved_labels = self.get_site_bond_resolved_labels()
+
+ for cation_name, icohp_data in bond_dict.items():
+ for atom_pair, bonding_data in dict_antibonding.items():
+ if (
+ namecation == atom_pair.split("-")[0]
+ and cation_name == atom_pair.split("-")[1]
+ ):
+ icohp_data["bonding"] = bonding_data["bonding"]
+ icohp_data["antibonding"] = bonding_data["antibonding"]
+ if self.orbital_resolved:
+ # get orb resolved data to be added
+ orb_resolved_bond_info = (
+ self._get_orbital_resolved_data(
+ nameion=namecation,
+ iion=ication,
+ labels=labels,
+ bond_resolved_labels=bond_resolved_labels,
+ type_pop=type_pop,
+ )
+ )
+ # match the dict key in bond_dict and get corresponding orbital data
+ for ion_atom_pair_orb in orb_resolved_bond_info:
+ orb_data_atom_pair = ion_atom_pair_orb.split(": ")[
+ -1
+ ]
+ atom_pair_here = atom_pair.split("-")
+ atom_pair_here.sort()
+ if (
+ orb_data_atom_pair == "-".join(atom_pair_here)
+ and (namecation + str(ication + 1) + ":")
+ in ion_atom_pair_orb
+ ):
+ icohp_data[
+ "orbital_data"
+ ] = orb_resolved_bond_info[ion_atom_pair_orb]
+
+ orb_data_stats = self._get_bond_resolved_data_stats(
+ orb_resolved_bond_data=orb_resolved_bond_info[
+ ion_atom_pair_orb
+ ],
+ )
+
+ icohp_data["orbital_data"].update(
+ orb_data_stats
+ )
+
+ site_dict[ication] = {
+ "env": ce,
+ "bonds": bond_dict,
+ "ion": namecation,
+ "charge": charge_list[ication],
+ "relevant_bonds": cation_anion_infos[3],
+ }
+ elif self.which_bonds == "all":
+ for iion, ce, bond_infos, labels, cohps in zip(
+ self.seq_ineq_ions,
+ self.seq_coord_ions,
+ self.seq_infos_bonds,
+ self.seq_labels_cohps,
+ self.seq_cohps,
+ ):
+ nameion = str(self.structure[iion].specie)
+
+ # This will compute the mean strengths of ICOHPs
+ mean_icohps = self._get_strenghts_for_each_bond(
+ pairs=bond_infos[4], strengths=bond_infos[1], nameion=None
+ )
+ # pairs, strengths, nameion
+ # will collect if there are antibonding states present
+ antbdg = self._get_antibdg_states(cohps, labels, nameion=None)
+
+ dict_antibonding = (
+ self._integrate_antbdstates_below_efermi_for_set_cohps(
+ labels, cohps, nameion
+ )
+ )
+
+ bond_dict = self._get_bond_dict(
+ mean_icohps, antbdg, nameion=nameion, type_pop=type_pop
+ )
+ bond_resolved_labels = self.get_site_bond_resolved_labels()
+
+ for cation_name, icohp_data in bond_dict.items():
+ for atom_pair, bonding_data in dict_antibonding.items():
+ if (
+ nameion == atom_pair.split("-")[0]
+ and cation_name == atom_pair.split("-")[1]
+ ):
+ icohp_data["bonding"] = bonding_data["bonding"]
+ icohp_data["antibonding"] = bonding_data["antibonding"]
+ if self.orbital_resolved:
+ # get orb resolved data to be added
+ orb_resolved_bond_info = (
+ self._get_orbital_resolved_data(
+ nameion=nameion,
+ iion=iion,
+ labels=labels,
+ bond_resolved_labels=bond_resolved_labels,
+ type_pop=type_pop,
+ )
+ )
+ # match the dict key in bond_dict and get corresponding orbital data
+ for ion_atom_pair_orb in orb_resolved_bond_info:
+ orb_data_atom_pair = ion_atom_pair_orb.split(": ")[
+ -1
+ ]
+ atom_pair_here = atom_pair.split("-")
+ atom_pair_here.sort()
+ if (
+ orb_data_atom_pair == "-".join(atom_pair_here)
+ and (nameion + str(iion + 1) + ":")
+ in ion_atom_pair_orb
+ ):
+ icohp_data[
+ "orbital_data"
+ ] = orb_resolved_bond_info[ion_atom_pair_orb]
+
+ orb_data_stats = self._get_bond_resolved_data_stats(
+ orb_resolved_bond_data=orb_resolved_bond_info[
+ ion_atom_pair_orb
+ ],
+ )
+
+ icohp_data["orbital_data"].update(
+ orb_data_stats
+ )
+
+ site_dict[iion] = {
+ "env": ce,
+ "bonds": bond_dict,
+ "ion": nameion,
+ "charge": charge_list[iion],
+ "relevant_bonds": bond_infos[3],
+ }
+
+ if self.path_to_madelung is None:
+ if self.which_bonds == "cation-anion":
+ # This sets the dictionary including the most important information on the compound
+ self.condensed_bonding_analysis = {
+ "formula": formula,
+ "max_considered_bond_length": max_bond_lengths,
+ f"limit_i{type_pop.lower()}": limit_icohps,
+ "number_of_considered_ions": number_considered_ions,
+ "sites": site_dict,
+ "type_charges": self.type_charge,
+ }
+ elif self.which_bonds == "all":
+ self.condensed_bonding_analysis = {
+ "formula": formula,
+ "max_considered_bond_length": max_bond_lengths,
+ f"limit_i{type_pop.lower()}": limit_icohps,
+ "number_of_considered_ions": number_considered_ions,
+ "sites": site_dict,
+ "type_charges": self.type_charge,
+ }
+ else:
+ from pymatgen.io.lobster import MadelungEnergies
+
+ madelung = MadelungEnergies(self.path_to_madelung)
+ if self.type_charge == "Mulliken":
+ madelung_energy = madelung.madelungenergies_Mulliken
+ elif self.type_charge == "Löwdin":
+ madelung_energy = madelung.madelungenergies_Loewdin
+ # This sets the dictionary including the most important information on the compound
+ if self.which_bonds == "cation-anion":
+ self.condensed_bonding_analysis = {
+ "formula": formula,
+ "max_considered_bond_length": max_bond_lengths,
+ f"limit_i{type_pop.lower()}": limit_icohps,
+ "number_of_considered_ions": number_considered_ions,
+ "sites": site_dict,
+ "type_charges": self.type_charge,
+ "madelung_energy": madelung_energy,
+ }
+ elif self.which_bonds == "all":
+ self.condensed_bonding_analysis = {
+ "formula": formula,
+ "max_considered_bond_length": max_bond_lengths,
+ f"limit_i{type_pop.lower()}": limit_icohps,
+ "number_of_considered_ions": number_considered_ions,
+ "sites": site_dict,
+ "type_charges": self.type_charge,
+ "madelung_energy": madelung_energy,
+ }
+
+[docs] def set_summary_dicts(self):
+ """
+ Set summary dict that can be used for correlations.
+
+ bond_dict that includes information on each bond
+
+ "has_antbd" tells if there are antbonding states
+ "ICOHP_mean" shows the mean of all ICOHPs in EV
+
+ {'Yb-Sb': { 'has_antbdg': False, 'ICOHP_mean': -1.7448},
+ 'Mn-Sb': { 'has_antbdg': True, 'ICOHP_mean': -1.525}}
+
+ a cation dict that includes all different coordination environments and counts for them
+ {'Na': {'T:4': 4, 'A:2': 4}, 'Si': {'T:6': 4, 'PP:6': 4}}
+
+ Returns:
+ None
+
+ """
+ relevant_ion_ids = [
+ isite for isite in self.list_equivalent_sites if isite in self.seq_ineq_ions
+ ]
+ # set population type
+ type_pop = self._get_pop_type()
+
+ final_dict_bonds = {}
+ for key in relevant_ion_ids:
+ item = self.condensed_bonding_analysis["sites"][key]
+ for type, properties in item["bonds"].items():
+ label_list = [item["ion"], str(type)]
+ new_label = sorted(label_list.copy())
+ label = str(new_label[0]) + "-" + str(new_label[1])
+
+ if label not in final_dict_bonds:
+ final_dict_bonds[label] = {
+ "number_of_bonds": int(properties["number_of_bonds"]),
+ f"I{type_pop}_sum": float(properties[f"I{type_pop}_sum"]),
+ "has_antbdg": properties["has_antibdg_states_below_Efermi"],
+ }
+ else:
+ final_dict_bonds[label]["number_of_bonds"] += int(
+ properties["number_of_bonds"]
+ )
+ final_dict_bonds[label][f"I{type_pop}_sum"] += float(
+ properties[f"I{type_pop}_sum"]
+ )
+ final_dict_bonds[label]["has_antbdg"] = (
+ final_dict_bonds[label]["has_antbdg"]
+ or properties["has_antibdg_states_below_Efermi"]
+ )
+ self.final_dict_bonds = {}
+ for key, item in final_dict_bonds.items():
+ self.final_dict_bonds[key] = {}
+ self.final_dict_bonds[key][f"I{type_pop}_mean"] = item[
+ f"I{type_pop}_sum"
+ ] / (item["number_of_bonds"])
+ self.final_dict_bonds[key]["has_antbdg"] = item["has_antbdg"]
+
+ # rework, add all environments!
+ final_dict_ions = {}
+ for key in relevant_ion_ids:
+ if (
+ self.condensed_bonding_analysis["sites"][key]["ion"]
+ not in final_dict_ions
+ ):
+ final_dict_ions[
+ self.condensed_bonding_analysis["sites"][key]["ion"]
+ ] = [self.condensed_bonding_analysis["sites"][key]["env"]]
+ else:
+ final_dict_ions[
+ self.condensed_bonding_analysis["sites"][key]["ion"]
+ ].append(self.condensed_bonding_analysis["sites"][key]["env"])
+
+ self.final_dict_ions = {}
+ for key, item in final_dict_ions.items():
+ self.final_dict_ions[key] = dict(Counter(item))
+
+[docs] @staticmethod
+ def get_lobster_calc_quality_summary(
+ path_to_poscar: str,
+ path_to_lobsterout: str,
+ path_to_lobsterin: str,
+ path_to_potcar: str | None = None,
+ potcar_symbols: list | None = None,
+ path_to_charge: str | None = None,
+ path_to_bandoverlaps: str | None = None,
+ path_to_doscar: str | None = None,
+ path_to_vasprun: str | None = None,
+ dos_comparison: bool = False,
+ e_range: list = [-5, 0],
+ n_bins: int | None = None,
+ bva_comp: bool = False,
+ ) -> dict:
+ """
+ Analyze LOBSTER calculation quality.
+
+ Args:
+ path_to_poscar: path to structure file
+ path_to_lobsterout: path to lobsterout file
+ path_to_lobsterin: path to lobsterin file
+ path_to_potcar: path to VASP potcar file
+ potcar_symbols: list of potcar symbols from postcar file (can be used if no potcar available)
+ path_to_charge: path to CHARGE.lobster file
+ path_to_bandoverlaps: path to bandOverlaps.lobster file
+ path_to_doscar: path to DOSCAR.lobster or DOSCAR.LSO.lobster file
+ path_to_vasprun: path to vasprun.xml file
+ dos_comparison: will compare DOS from VASP and LOBSTER and return tanimoto index
+ e_range: energy range for DOS comparisons
+ n_bins: number of bins to discretize DOS for comparisons
+ bva_comp: Compares LOBSTER charge signs with Bond valence charge signs
+
+ Returns:
+ A dict of summary of LOBSTER calculation quality by analyzing basis set used,
+ charge spilling from lobsterout/ PDOS comparisons of VASP and LOBSTER /
+ BVA charge comparisons
+
+ """
+ quality_dict = {}
+
+ if path_to_potcar and not potcar_symbols:
+ potcar_names = Lobsterin._get_potcar_symbols(POTCAR_input=path_to_potcar)
+ elif not path_to_potcar and potcar_symbols:
+ potcar_names = potcar_symbols
+ else:
+ raise ValueError(
+ "Please provide either path_to_potcar or list of "
+ "potcar_symbols used for the calculations"
+ )
+
+ struct = Structure.from_file(path_to_poscar)
+
+ ref_bases = Lobsterin.get_all_possible_basis_functions(
+ structure=struct, potcar_symbols=potcar_names
+ )
+
+ lobs_in = Lobsterin.from_file(path_to_lobsterin)
+ calc_basis = []
+ for basis in lobs_in["basisfunctions"]:
+ basis_sep = basis.split()[1:]
+ basis_comb = " ".join(basis_sep)
+ calc_basis.append(basis_comb)
+
+ if calc_basis == list(ref_bases[0].values()):
+ quality_dict["minimal_basis"] = True # type: ignore
+ else:
+ quality_dict["minimal_basis"] = False # type: ignore
+ warnings.warn(
+ "Consider rerunning the calc with the minimum basis as well. Choosing is "
+ "larger basis set is recommended if you see a significant improvement of "
+ "the charge spilling and material has non-zero band gap."
+ )
+
+ lob_out = Lobsterout(path_to_lobsterout)
+
+ quality_dict["charge_spilling"] = {
+ "abs_charge_spilling": round((sum(lob_out.charge_spilling) / 2) * 100, 4),
+ "abs_total_spilling": round((sum(lob_out.total_spilling) / 2) * 100, 4),
+ } # type: ignore
+
+ if path_to_bandoverlaps is not None:
+ if Path(path_to_bandoverlaps).exists(): # type: ignore
+ band_overlaps = Bandoverlaps(filename=path_to_bandoverlaps)
+ for line in lob_out.warning_lines:
+ if "k-points could not be orthonormalized" in line:
+ total_kpoints = int(line.split(" ")[2])
+
+ # store actual number of devations above pymatgen default limit of 0.1
+ dev_val = []
+ for dev in band_overlaps.max_deviation:
+ if dev > 0.1:
+ dev_val.append(dev)
+
+ quality_dict["band_overlaps"] = { # type: ignore
+ "file_exists": True,
+ "limit_maxDeviation": 0.1,
+ "has_good_quality_maxDeviation": band_overlaps.has_good_quality_maxDeviation(
+ limit_maxDeviation=0.1
+ ),
+ "max_deviation": round(max(band_overlaps.max_deviation), 4),
+ "percent_kpoints_abv_limit": round(
+ (len(dev_val) / total_kpoints) * 100, 4
+ ),
+ }
+
+ else:
+ quality_dict["band_overlaps"] = {
+ "file_exists": False,
+ "limit_maxDeviation": None,
+ "has_good_quality_maxDeviation": True,
+ "max_deviation": None,
+ "percent_kpoints_abv_limit": None,
+ } # type: ignore
+
+ if bva_comp:
+ try:
+ bond_valence = BVAnalyzer()
+
+ bva_oxi = []
+ lobs_charge = Charge(filename=path_to_charge)
+ for i in bond_valence.get_valences(structure=struct):
+ if i >= 0:
+ bva_oxi.append("POS")
+ else:
+ bva_oxi.append("NEG")
+
+ mull_oxi = []
+ for i in lobs_charge.Mulliken:
+ if i >= 0:
+ mull_oxi.append("POS")
+ else:
+ mull_oxi.append("NEG")
+
+ loew_oxi = []
+ for i in lobs_charge.Loewdin:
+ if i >= 0:
+ loew_oxi.append("POS")
+ else:
+ loew_oxi.append("NEG")
+
+ quality_dict["charges"] = {} # type: ignore
+ if mull_oxi == bva_oxi:
+ quality_dict["charges"]["BVA_Mulliken_agree"] = True # type: ignore
+ else:
+ quality_dict["charges"]["BVA_Mulliken_agree"] = False # type: ignore
+
+ if mull_oxi == bva_oxi:
+ quality_dict["charges"]["BVA_Loewdin_agree"] = True # type: ignore
+ else:
+ quality_dict["charges"]["BVA_Loewdin_agree"] = False # type: ignore
+ except ValueError:
+ quality_dict["charges"] = {} # type: ignore
+ warnings.warn(
+ "Oxidation states from BVA analyzer cannot be determined. "
+ "Thus BVA charge comparison will be skipped"
+ )
+ if dos_comparison:
+ if "LSO" not in str(path_to_doscar).split("."):
+ warnings.warn(
+ "Consider using DOSCAR.LSO.lobster, as non LSO DOS from LOBSTER can have "
+ "negative DOS values"
+ )
+ doscar_lobster = Doscar(
+ doscar=path_to_doscar,
+ structure_file=path_to_poscar,
+ )
+
+ dos_lobster = doscar_lobster.completedos
+
+ vasprun = Vasprun(path_to_vasprun)
+ dos_vasp = vasprun.complete_dos
+
+ quality_dict["dos_comparisons"] = {} # type: ignore
+
+ for orb in dos_lobster.get_spd_dos():
+ if e_range[0] >= min(dos_vasp.energies) and e_range[0] >= min(
+ dos_lobster.energies
+ ):
+ min_e = e_range[0]
+ else:
+ warnings.warn(
+ "Minimum energy range requested for DOS comparisons is not available "
+ "in VASP or LOBSTER calculation. Thus, setting min_e to -5 eV"
+ )
+ min_e = -5
+
+ if e_range[-1] <= max(dos_vasp.energies) and e_range[-1] <= max(
+ dos_lobster.energies
+ ):
+ max_e = e_range[-1]
+ else:
+ warnings.warn(
+ "Maximum energy range requested for DOS comparisons is not available "
+ "in VASP or LOBSTER calculation. Thus, setting max_e to 0 eV"
+ )
+ max_e = 0
+
+ if (
+ np.diff(dos_vasp.energies)[0] >= 0.1
+ or np.diff(dos_lobster.energies)[0] >= 0.1
+ ):
+ warnings.warn(
+ "Input DOS files have very few points in the energy interval and thus "
+ "comparisons will not be reliable. Please rerun the calculations with "
+ "higher number of DOS points. Set NEDOS and COHPSteps tags to >= 2000 in VASP and LOBSTER "
+ "calculations, respectively."
+ )
+
+ if not n_bins:
+ n_bins = 56
+
+ fp_lobster_orb = dos_lobster.get_dos_fp(
+ min_e=min_e,
+ max_e=max_e,
+ n_bins=n_bins,
+ normalize=True,
+ type=orb.name,
+ )
+ fp_vasp_orb = dos_vasp.get_dos_fp(
+ min_e=min_e,
+ max_e=max_e,
+ n_bins=n_bins,
+ normalize=True,
+ type=orb.name,
+ )
+
+ tani_orb = round(
+ dos_vasp.get_dos_fp_similarity(
+ fp_lobster_orb, fp_vasp_orb, tanimoto=True
+ ),
+ 4,
+ )
+ quality_dict["dos_comparisons"][
+ f"tanimoto_orb_{orb.name}"
+ ] = tani_orb # type: ignore
+
+ fp_lobster = dos_lobster.get_dos_fp(
+ min_e=min_e,
+ max_e=max_e,
+ n_bins=n_bins,
+ normalize=True,
+ type="summed_pdos",
+ )
+ fp_vasp = dos_vasp.get_dos_fp(
+ min_e=min_e,
+ max_e=max_e,
+ n_bins=n_bins,
+ normalize=True,
+ type="summed_pdos",
+ )
+
+ tanimoto_summed = round(
+ dos_vasp.get_dos_fp_similarity(fp_lobster, fp_vasp, tanimoto=True), 4
+ )
+ quality_dict["dos_comparisons"]["tanimoto_summed"] = tanimoto_summed # type: ignore
+ quality_dict["dos_comparisons"]["e_range"] = [min_e, max_e] # type: ignore
+ quality_dict["dos_comparisons"]["n_bins"] = n_bins # type: ignore
+
+ return quality_dict
+
+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""This module defines classes to describe the COHPs automatically."""
+from __future__ import annotations
+
+from pathlib import Path
+
+from lobsterpy.plotting import InteractiveCohpPlotter, PlainCohpPlotter
+
+
+[docs]class Description:
+ """
+ Base class that will write generate a text description for all relevant bonds.
+
+ It analyses all relevant coordination environments in the system based on electronic structure theory.
+
+ """
+
+ def __init__(self, analysis_object):
+ """
+ Generate a text description for all relevant bonds.
+
+ Args:
+ analysis_object: Analysis object from lobsterpy.analysis
+ """
+ self.analysis_object = analysis_object
+ self.set_description()
+
+[docs] def set_description(self):
+ """
+ Set the descriptions of the structures using the cation names, starting with numbers at 1.
+
+ Uses the cation names from the lobster files.
+
+ Returns:
+ None
+
+ """
+ self.condensed_bonding_analysis = (
+ self.analysis_object.condensed_bonding_analysis
+ )
+ # set type of population analyzed
+ type_pop = self.analysis_object._get_pop_type()
+ # set units for populations
+ units = " eV" if type_pop == "COHP" else ""
+ if self.analysis_object.which_bonds == "cation-anion":
+ relevant_cations = ", ".join(
+ [
+ str(site.specie) + str(isite + 1)
+ for isite, site in enumerate(self.analysis_object.structure)
+ if isite in self.analysis_object.seq_ineq_ions
+ ]
+ )
+ self.text = []
+ self.text.append(
+ "The compound "
+ + str(self.condensed_bonding_analysis["formula"])
+ + " has "
+ + str(self.condensed_bonding_analysis["number_of_considered_ions"])
+ + " symmetry-independent cation(s) with relevant cation-anion interactions: "
+ + str(relevant_cations)
+ + "."
+ )
+
+ for key, item in self.condensed_bonding_analysis["sites"].items():
+ # It has 3 Ta-N (mean ICOHP: -4.78 eV, antibonding interactions below EFermi),
+ bond_info = []
+ orb_info = []
+ for type, properties in item["bonds"].items():
+ if not properties["has_antibdg_states_below_Efermi"]:
+ bond_info.append(
+ str(properties["number_of_bonds"])
+ + " "
+ + item["ion"]
+ + "-"
+ + str(type)
+ + f" (mean I{type_pop}: "
+ ""
+ + properties[f"I{type_pop}_mean"]
+ + f"{units}, 0.0 percent antibonding interaction below EFermi)"
+ )
+ if self.analysis_object.orbital_resolved:
+ text_orbital = (
+ self._generate_orbital_resolved_analysis_text(
+ orbital_resolved_data=properties,
+ type_pop=type_pop,
+ atom_name=str(type),
+ ion=item["ion"],
+ )
+ )
+ orb_info.extend(text_orbital)
+ else:
+ bond_info.append(
+ str(properties["number_of_bonds"])
+ + " "
+ + item["ion"]
+ + "-"
+ + str(type)
+ + f" (mean I{type_pop}: "
+ ""
+ + properties[f"I{type_pop}_mean"]
+ + f"{units}, "
+ + str(round(properties["antibonding"]["perc"] * 100, 3))
+ + " percent antibonding interaction below EFermi)"
+ )
+ if self.analysis_object.orbital_resolved:
+ text_orbital = (
+ self._generate_orbital_resolved_analysis_text(
+ orbital_resolved_data=properties,
+ type_pop=type_pop,
+ atom_name=str(type),
+ ion=item["ion"],
+ )
+ )
+ orb_info.extend(text_orbital)
+
+ bonds = (
+ ",".join(bond_info[0:-1]) + ", and " + bond_info[-1]
+ if len(bond_info) > 1
+ else bond_info[0]
+ )
+
+ if len(orb_info) > 1:
+ orb_bonds = "".join(orb_info).replace(".In", ". In")
+ else:
+ orb_bonds = orb_info[0] if orb_info else ""
+ if item["env"] == "O:6":
+ self.text.append(
+ str(item["ion"])
+ + str(key + 1)
+ + " has an "
+ + str(self._coordination_environment_to_text(item["env"]))
+ + " coordination environment. It has "
+ + str(bonds)
+ + " bonds."
+ )
+ if orb_bonds:
+ self.text.append(orb_bonds)
+ else:
+ self.text.append(
+ str(item["ion"])
+ + str(key + 1)
+ + " has a "
+ + str(self._coordination_environment_to_text(item["env"]))
+ + " coordination environment. It has "
+ + str(bonds)
+ + " bonds."
+ )
+ if orb_bonds:
+ self.text.append(orb_bonds)
+
+ elif self.analysis_object.which_bonds == "all":
+ relevant_ions = ", ".join(
+ [
+ str(site.specie) + str(isite + 1)
+ for isite, site in enumerate(self.analysis_object.structure)
+ if isite in self.analysis_object.seq_ineq_ions
+ ]
+ )
+ self.text = []
+ self.text.append(
+ "The compound "
+ + str(self.condensed_bonding_analysis["formula"])
+ + " has "
+ + str(self.condensed_bonding_analysis["number_of_considered_ions"])
+ + " symmetry-independent atoms(s) with relevant bonds: "
+ + str(relevant_ions)
+ + "."
+ )
+
+ for key, item in self.condensed_bonding_analysis["sites"].items():
+ # It has 3 Ta-N (mean ICOHP: -4.78 eV, antibonding interactions below EFermi),
+ bond_info = []
+ orb_info = []
+ for type, properties in item["bonds"].items():
+ if not properties["has_antibdg_states_below_Efermi"]:
+ bond_info.append(
+ str(properties["number_of_bonds"])
+ + " "
+ + item["ion"]
+ + "-"
+ + str(type)
+ + f" (mean I{type_pop}: "
+ ""
+ + properties[f"I{type_pop}_mean"]
+ + f"{units}, 0.0 percent antibonding interaction below EFermi)"
+ )
+ if self.analysis_object.orbital_resolved:
+ text_orbital = (
+ self._generate_orbital_resolved_analysis_text(
+ orbital_resolved_data=properties,
+ type_pop=type_pop,
+ atom_name=str(type),
+ ion=item["ion"],
+ )
+ )
+ orb_info.extend(text_orbital)
+ else:
+ bond_info.append(
+ str(properties["number_of_bonds"])
+ + " "
+ + item["ion"]
+ + "-"
+ + str(type)
+ + f" (mean I{type_pop}: "
+ ""
+ + properties[f"I{type_pop}_mean"]
+ + f"{units}, "
+ + str(round(properties["antibonding"]["perc"] * 100, 3))
+ + " percent antibonding interaction below EFermi)"
+ )
+
+ if self.analysis_object.orbital_resolved:
+ text_orbital = (
+ self._generate_orbital_resolved_analysis_text(
+ orbital_resolved_data=properties,
+ type_pop=type_pop,
+ atom_name=str(type),
+ ion=item["ion"],
+ )
+ )
+ orb_info.extend(text_orbital)
+
+ if len(bond_info) > 1:
+ bonds = ",".join(bond_info[0:-1]) + ", and " + bond_info[-1]
+ else:
+ bonds = bond_info[0] if bond_info else 0
+ if len(orb_info) > 1:
+ orb_bonds = "".join(orb_info).replace(".In", ". In")
+ else:
+ orb_bonds = orb_info[0] if orb_info else ""
+ if item["env"] == "O:6":
+ self.text.append(
+ str(item["ion"])
+ + str(key + 1)
+ + " has an "
+ + str(self._coordination_environment_to_text(item["env"]))
+ + " coordination environment. It has "
+ + str(bonds)
+ + " bonds."
+ )
+ if orb_bonds:
+ self.text.append(orb_bonds)
+ else:
+ self.text.append(
+ str(item["ion"])
+ + str(key + 1)
+ + " has a "
+ + str(self._coordination_environment_to_text(item["env"]))
+ + " coordination environment. It has "
+ + str(bonds)
+ + " bonds."
+ )
+ if orb_bonds:
+ self.text.append(orb_bonds)
+
+ if "madelung_energy" in self.analysis_object.condensed_bonding_analysis:
+ self.text.append(
+ "The Madelung energy of this crystal structure per unit cell is: "
+ + str(
+ self.analysis_object.condensed_bonding_analysis["madelung_energy"]
+ )
+ + " eV."
+ )
+
+[docs] def _generate_orbital_resolved_analysis_text(
+ self,
+ orbital_resolved_data: dict,
+ ion: str,
+ atom_name: str,
+ type_pop: str,
+ ):
+ """
+ Generate text from orbital-resolved analysis data of the most relevant COHP, COOP, or COBI.
+
+ Args:
+ orbital_resolved_data : dict of orbital data from condensed bonding analysis object
+ ion: name of ion at the site
+ atom_name: name of atomic speice to which ion is bonded
+ type_pop: population type analysed could be "COHP" or "COOP" or "COBI"
+
+ Returns:
+ A python list with text describing the orbital which contributes
+ the most to the bonding and antibonding in the bond at site
+ """
+ orb_info = []
+ if orbital_resolved_data["orbital_data"]["orbital_summary_stats"]:
+ orb_names = []
+ orb_contri = []
+ # get atom-pair list with ion placed first
+ atom_pair = self.analysis_object._sort_name([ion, atom_name], nameion=ion)
+ if (
+ "max_bonding_contribution"
+ in orbital_resolved_data["orbital_data"]["orbital_summary_stats"]
+ ):
+ for orb, data in orbital_resolved_data["orbital_data"][
+ "orbital_summary_stats"
+ ]["max_bonding_contribution"].items():
+ atom_pair_with_orb_name = (
+ self.analysis_object._sort_orbital_atom_pair(
+ atom_pair=atom_pair,
+ complete_cohp=self.analysis_object.chemenv.completecohp,
+ label=orbital_resolved_data["orbital_data"][
+ "relevant_bonds"
+ ][0],
+ orb_pair=orb,
+ )
+ )
+ orb_names.append("-".join(atom_pair_with_orb_name))
+ orb_contri.append(
+ str(
+ round(
+ data * 100,
+ 3,
+ )
+ )
+ )
+ orb_names_anti = []
+ orb_antibonding = []
+ if (
+ "max_antibonding_contribution"
+ in orbital_resolved_data["orbital_data"]["orbital_summary_stats"]
+ ):
+ for orb, data in orbital_resolved_data["orbital_data"][
+ "orbital_summary_stats"
+ ]["max_antibonding_contribution"].items():
+ atom_pair_with_orb_name = (
+ self.analysis_object._sort_orbital_atom_pair(
+ atom_pair=atom_pair,
+ complete_cohp=self.analysis_object.chemenv.completecohp,
+ label=orbital_resolved_data["orbital_data"][
+ "relevant_bonds"
+ ][0],
+ orb_pair=orb,
+ )
+ )
+ orb_names_anti.append("-".join(atom_pair_with_orb_name))
+ orb_antibonding.append(
+ str(
+ round(
+ data * 100,
+ 3,
+ )
+ )
+ )
+ if len(orb_contri) > 1:
+ orb_name_contri = ""
+ for inx, name in enumerate(orb_names):
+ if len(orb_contri) == 2 and inx + 1 != len(orb_contri):
+ orb_name_contri += f"{name} "
+ elif 2 < len(orb_contri) != inx + 1:
+ orb_name_contri += f"{name}, "
+ else:
+ orb_name_contri += f"and {name}"
+
+ orb_name_contri += " orbitals, contributing "
+ for inx, contribution in enumerate(orb_contri):
+ if len(orb_contri) == 2 and inx + 1 != len(orb_contri):
+ orb_name_contri += f"{contribution} "
+ elif 2 < len(orb_contri) != inx + 1:
+ orb_name_contri += f"{contribution}, "
+ else:
+ orb_name_contri += f"and {contribution} percent, respectively"
+ num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"])
+ bonds = "bonds" if num_bonds > 1 else "bond"
+ orb_info.append(
+ f"In the {num_bonds} "
+ + "-".join(atom_pair)
+ + f" {bonds}, relative to the summed I{type_pop}s, "
+ + "the maximum bonding contribution is from "
+ + orb_name_contri
+ )
+ elif not orb_contri:
+ num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"])
+ bonds = "bonds" if num_bonds > 1 else "bond"
+ orb_info.append(
+ f"In the {num_bonds} "
+ + "-".join(atom_pair)
+ + f" {bonds}, relative to the summed I{type_pop}s, "
+ + f"no orbital has a bonding contribution greater than "
+ f"{self.analysis_object.orbital_cutoff*100} percent"
+ )
+ else:
+ num_bonds = len(orbital_resolved_data["orbital_data"]["relevant_bonds"])
+ bonds = "bonds" if num_bonds > 1 else "bond"
+ orb_info.append(
+ f"In the {num_bonds} "
+ + "-".join(atom_pair)
+ + f" {bonds}, relative to the summed I{type_pop}s, "
+ + "the maximum bonding contribution is from the "
+ + f"{orb_names[0]}"
+ + f" orbital, contributing {orb_contri[0]} percent"
+ )
+
+ if len(orb_antibonding) > 1:
+ orb_anti = ""
+ for inx, name in enumerate(orb_names_anti):
+ if len(orb_names_anti) == 2 and inx + 1 != len(orb_names_anti):
+ orb_anti += f"{name} "
+ elif 2 < len(orb_antibonding) != inx + 1:
+ orb_anti += f"{name}, "
+ else:
+ orb_anti += f"and {name}"
+
+ orb_anti += " orbitals, contributing "
+ for inx, contribution in enumerate(orb_antibonding):
+ if len(orb_names_anti) == 2 and inx + 1 != len(orb_names_anti):
+ orb_anti += f"{contribution} "
+ elif 2 < len(orb_antibonding) != inx + 1:
+ orb_anti += f"{contribution}, "
+ else:
+ orb_anti += f"and {contribution} percent, respectively."
+ orb_info.append(
+ f", whereas the maximum antibonding contribution is from {orb_anti}"
+ )
+ elif not orb_antibonding:
+ orb_info.append(
+ ", whereas no significant antibonding contribution is found in this bond."
+ )
+ else:
+ orb_info.append(
+ f", whereas the maximum antibonding contribution is from the "
+ f"{orb_names_anti[0]} orbital, contributing {orb_antibonding[0]} percent."
+ )
+ else:
+ # get atom-pair list with ion placed first
+ atom_pair = self.analysis_object._sort_name([ion, atom_name], nameion=ion)
+ percentage_cutoff = round(self.analysis_object.orbital_cutoff * 100, 2)
+ orb_info.append(
+ f"No individual orbital interactions detected above {percentage_cutoff} percent"
+ f" with summed I{type_pop} as reference for the "
+ + "-".join(atom_pair)
+ + " bond."
+ )
+
+ return orb_info
+
+[docs] def plot_cohps(
+ self,
+ save=False,
+ filename=None,
+ ylim=[-4, 2],
+ xlim=None,
+ integrated=False,
+ title="",
+ sigma=None,
+ hide=False,
+ ):
+ """
+ Automatically generate plots of the most relevant COHPs, COOPs, or COBIs.
+
+ Args:
+ save (bool): will save the plot to a file
+ filename (str/Path): name of the file to save the plot.
+ ylim (list of float): energy scale that is shown in plot (eV)
+ xlim(list of float): energy range for COHPs in eV
+ integrated (bool): if True, integrated COHPs will be shown
+ sigma: Standard deviation of Gaussian broadening applied to
+ population data. If None, no broadening will be added.
+ title: sets the title of figure generated
+ hide (bool): if True, the plot will not be shown.
+
+ Returns:
+ A matplotlib object.
+
+ """
+ seq_cohps = self.analysis_object.seq_cohps
+ if self.analysis_object.which_bonds == "cation-anion":
+ seq_ineq_cations = self.analysis_object.seq_ineq_ions
+ elif self.analysis_object.which_bonds == "all":
+ seq_ineq_cations = self.analysis_object.seq_ineq_ions
+ seq_labels = self.analysis_object.seq_labels_cohps
+ structure = self.analysis_object.structure
+
+ for iplot, (ication, labels, cohps) in enumerate(
+ zip(seq_ineq_cations, seq_labels, seq_cohps)
+ ):
+ namecation = str(structure[ication].specie)
+
+ cp = PlainCohpPlotter(
+ are_coops=self.analysis_object.are_coops,
+ are_cobis=self.analysis_object.are_cobis,
+ )
+ for label, cohp in zip(labels, cohps):
+ if label is not None:
+ cp.add_cohp(namecation + str(ication + 1) + ": " + label, cohp)
+
+ plot = cp.get_plot(integrated=integrated, sigma=sigma)
+ plot.ylim(ylim)
+ if xlim is not None:
+ plot.xlim(xlim)
+
+ plot.title(title)
+ if save:
+ if len(seq_ineq_cations) > 1:
+ if isinstance(filename, str):
+ filename = Path(filename)
+ filename_new = (
+ filename.parent / f"{filename.stem}-{iplot}{filename.suffix}"
+ )
+ else:
+ filename_new = filename
+ plot.savefig(filename_new)
+ if not hide:
+ plot.show()
+ else:
+ if not save:
+ plot.close()
+
+[docs] def plot_interactive_cohps(
+ self,
+ save_as_html=False,
+ filename=None,
+ ylim=None,
+ xlim=None,
+ integrated=False,
+ title="",
+ sigma=None,
+ label_resolved=False,
+ orbital_resolved=False,
+ hide=False,
+ ):
+ """
+ Automatically generate interactive plots of the most relevant COHPs, COBIs or COOPs.
+
+ Args:
+ save_as_html (bool): will save the plot to a html file
+ filename (str/Path): name of the file to save the plot.
+ ylim (list of float): energy scale that is shown in plot (eV)
+ xlim (list of float): energy range for COHPs in eV
+ integrated (bool): if True, integrated COHPs will be shown
+ sigma: Standard deviation of Gaussian broadening applied to
+ population data. If None, no broadening will be added.
+ label_resolved: if true, relevant cohp curves will be further resolved based on band labels
+ orbital_resolved: if true, relevant orbital interactions in cohp curves will be added to figure
+ title : Title of the interactive plot
+ hide (bool): if True, the plot will not be shown.
+
+ Returns:
+ A plotly.graph_objects.Figure object.
+ """
+ cba_cohp_plot_data = {} # Initialize dict to store plot data
+ set_cohps = self.analysis_object.seq_cohps
+ set_labels_cohps = self.analysis_object.seq_labels_cohps
+ set_inequivalent_cations = self.analysis_object.seq_ineq_ions
+ structure = self.analysis_object.structure
+
+ for _iplot, (ication, labels, cohps) in enumerate(
+ zip(set_inequivalent_cations, set_labels_cohps, set_cohps)
+ ):
+ label_str = f"{structure[ication].specie!s}{ication + 1!s}: "
+ for label, cohp in zip(labels, cohps):
+ if label is not None:
+ cba_cohp_plot_data[label_str + label] = cohp
+
+ ip = InteractiveCohpPlotter(
+ are_coops=self.analysis_object.are_coops,
+ are_cobis=self.analysis_object.are_cobis,
+ )
+ if label_resolved or orbital_resolved:
+ ip.add_all_relevant_cohps(
+ analyse=self.analysis_object,
+ label_resolved=label_resolved,
+ orbital_resolved=orbital_resolved,
+ )
+ else:
+ ip.add_cohps_from_plot_data(plot_data_dict=cba_cohp_plot_data)
+
+ plot = ip.get_plot(integrated=integrated, xlim=xlim, ylim=ylim, sigma=sigma)
+
+ plot.update_layout(title_text=title)
+ if save_as_html:
+ plot.write_html(filename, include_mathjax="cdn")
+ if not hide:
+ return plot.show()
+
+ return plot
+
+[docs] @staticmethod
+ def _coordination_environment_to_text(ce):
+ """
+ Convert a coordination environment string into a text description of the environment.
+
+ Args:
+ ce (str): output from ChemEnv package (e.g., "O:6")
+
+ Returns:
+ A text description of coordination environment
+ """
+ if ce == "S:1":
+ return "single (CN=1)"
+ if ce == "L:2":
+ return "linear (CN=2)"
+ if ce == "A:2":
+ return "angular (CN=2)"
+ if ce == "TL:3":
+ return "trigonal planar (CN=3)"
+ if ce == "TY:3":
+ return "triangular non-coplanar (CN=3)"
+ if ce == "TS:3":
+ return "t-shaped (CN=3)"
+ if ce == "T:4":
+ return "tetrahedral (CN=4)"
+ if ce == "S:4":
+ return "square planar (CN=4)"
+ if ce == "SY:4":
+ return "square non-coplanar (CN=4)"
+ if ce == "SS:4":
+ return "see-saw like (CN=4)"
+ if ce == "PP:5":
+ return "pentagonal (CN=5)"
+ if ce == "S:5":
+ return "square pyramidal (CN=5)"
+ if ce == "T:5":
+ return "trigonal bipyramidal (CN=5)"
+ if ce == "O:6":
+ return "octahedral (CN=6)"
+ if ce == "T:6":
+ return "trigonal prismatic (CN=6)"
+ if ce == "PP:6":
+ return "pentagonal pyramidal (CN=6)"
+ if ce == "PB:7":
+ return "pentagonal bipyramidal (CN=7)"
+ if ce == "ST:7":
+ return "square-face capped trigonal prismatic (CN=7)"
+ if ce == "ET:7":
+ return "end-trigonal-face capped trigonal prismatic (CN=7)"
+ if ce == "FO:7":
+ return "face-capped octahedron (CN=7)"
+ if ce == "C:8":
+ return "cubic (CN=8)"
+ if ce == "SA:8":
+ return "square antiprismatic (CN=8)"
+ if ce == "SBT:8":
+ return "square-face bicapped trigonal prismatic (CN=8)"
+ if ce == "TBT:8":
+ return "triangular-face bicapped trigonal prismatic (CN=8)"
+ if ce == "DD:8":
+ return "dodecahedronal (with triangular faces) (CN=8)"
+ if ce == "DDPN:8":
+ return (
+ "dodecahedronal (with triangular faces - p2345 plane normalized) (CN=8)"
+ )
+ if ce == "HB:8":
+ return "hexagonal bipyramidal (CN=8)"
+ if ce == "BO_1:8":
+ return "bicapped octahedral (opposed cap faces) (CN=8)"
+ if ce == "BO_2:8":
+ return "bicapped octahedral (cap faces with one atom in common) (CN=8)"
+ if ce == "BO_3:8":
+ return "bicapped octahedral (cap faces with one edge in common) (CN=8)"
+ if ce == "TC:9":
+ return "triangular cupola (CN=9)"
+ if ce == "TT_1:9":
+ return "Tricapped triangular prismatic (three square - face caps) (CN=9)"
+ if ce == "TT_2:9":
+ return "Tricapped triangular prismatic (two square - face caps and one triangular - face cap) (CN=9)"
+ if ce == "TT_3:9":
+ return "Tricapped triangular prism (one square - face cap and two triangular - face caps) (CN=9)"
+ if ce == "HD:9":
+ return "Heptagonal dipyramidal (CN=9)"
+ if ce == "TI:9":
+ return "tridiminished icosohedral (CN=9)"
+ if ce == "SMA:9":
+ return "Square-face monocapped antiprism (CN=9)"
+ if ce == "SS:9":
+ return "Square-face capped square prismatic (CN=9)"
+ if ce == "TO_1:9":
+ return "Tricapped octahedral (all 3 cap faces share one atom) (CN=9)"
+ if ce == "TO_2:9":
+ return "Tricapped octahedral (cap faces are aligned) (CN=9)"
+ if ce == "TO_3:9":
+ return "Tricapped octahedron (all 3 cap faces are sharing one edge of a face) (CN=9)"
+ if ce == "PP:10":
+ return "Pentagonal prismatic (CN=10)"
+ if ce == "PA:10":
+ return "Pentagonal antiprismatic (CN=10)"
+ if ce == "SBSA:10":
+ return "Square-face bicapped square antiprismatic (CN=10)"
+ if ce == "MI:10":
+ return "Metabidiminished icosahedral (CN=10)"
+ if ce == "S:10":
+ return "sphenocoronal (CN=10)"
+ if ce == "H:10":
+ return "Hexadecahedral (CN=10)"
+ if ce == "BS_1:10":
+ return "Bicapped square prismatic (opposite faces) (CN=10)"
+ if ce == "BS_1:10":
+ return "Bicapped square prismatic (opposite faces) (CN=10)"
+ if ce == "BS_2:10":
+ return "Bicapped square prism(adjacent faces) (CN=10)"
+ if ce == "TBSA:10":
+ return "Trigonal-face bicapped square antiprismatic (CN=10)"
+ if ce == "PCPA:11":
+ return "Pentagonal - face capped pentagonal antiprismatic (CN=11)"
+ if ce == "H:11":
+ return "Hendecahedral (CN=11)"
+ if ce == "SH:11":
+ return "Sphenoid hendecahedral (CN=11)"
+ if ce == "CO:11":
+ return "Cs - octahedral (CN=11)"
+ if ce == "DI:11":
+ return "Diminished icosahedral (CN=12)"
+ if ce == "I:12":
+ return "Icosahedral (CN=12)"
+ if ce == "PBP: 12":
+ return "Pentagonal - face bicapped pentagonal prismatic (CN=12)"
+ if ce == "TT:12":
+ return "Truncated tetrahedral (CN=12)"
+ if ce == "C:12":
+ return "Cuboctahedral (CN=12)"
+ if ce == "AC:12":
+ return "Anticuboctahedral (CN=12)"
+ if ce == "SC:12":
+ return "Square cupola (CN=12)"
+ if ce == "S:12":
+ return "Sphenomegacorona (CN=12)"
+ if ce == "HP:12":
+ return "Hexagonal prismatic (CN=12)"
+ if ce == "HA:12":
+ return "Hexagonal antiprismatic (CN=12)"
+ if ce == "SH:13":
+ return "Square-face capped hexagonal prismatic (CN=13)"
+ if ce == "1":
+ return "1-fold"
+ if ce == "2":
+ return "2-fold"
+ if ce == "3":
+ return "3-fold"
+ if ce == "4":
+ return "4-fold"
+ if ce == "5":
+ return "5-fold"
+ if ce == "6":
+ return "6-fold"
+ if ce == "7":
+ return "7-fold"
+ if ce == "8":
+ return "8-fold"
+ if ce == "9":
+ return "9-fold"
+ if ce == "10":
+ return "10-fold"
+ if ce == "11":
+ return "11-fold"
+ if ce == "12":
+ return "12-fold"
+ if ce == "13":
+ return "13-fold"
+ if ce == "14":
+ return "14-fold"
+ if ce == "15":
+ return "15-fold"
+ if ce == "16":
+ return "16-fold"
+ if ce == "17":
+ return "17-fold"
+ if ce == "18":
+ return "18-fold"
+ if ce == "19":
+ return "19-fold"
+ if ce == "20":
+ return "20-fold"
+ if ce == "21":
+ return "21-fold"
+ if ce == "22":
+ return "22-fold"
+ if ce == "23":
+ return "23-fold"
+ if ce == "24":
+ return "24-fold"
+ if ce == "25":
+ return "25-fold"
+ if ce == "26":
+ return "26-fold"
+ if ce == "27":
+ return "27-fold"
+ if ce == "28":
+ return "28-fold"
+ if ce == "29":
+ return "29-fold"
+ if ce == "30":
+ return "30-fold"
+ return ce
+
+[docs] def write_description(self):
+ """Print the description of the COHPs or COBIs or COOPs to the screen."""
+ for textpart in self.text:
+ print(textpart)
+
+[docs] @staticmethod
+ def get_calc_quality_description(quality_dict):
+ """
+ Generate a text description of the LOBSTER calculation quality.
+
+ Args:
+ quality_dict: python dictionary from lobsterpy.analysis.get_lobster_calc_quality_summary
+ """
+ text_des = []
+
+ for key, val in quality_dict.items():
+ if key == "minimal_basis":
+ if val:
+ text_des.append("The LOBSTER calculation used minimal basis.")
+ if not val:
+ text_des.append(
+ "Consider rerunning the calculation with the minimum basis as well. Choosing a "
+ "larger basis set is only recommended if you see a significant improvement of "
+ "the charge spilling."
+ )
+
+ elif key == "charge_spilling":
+ text_des.append(
+ "The absolute and total charge spilling for the calculation is {} and {} %, "
+ "respectively.".format(
+ quality_dict[key]["abs_charge_spilling"],
+ quality_dict[key]["abs_total_spilling"],
+ )
+ )
+ elif key == "band_overlaps":
+ if quality_dict[key]["file_exists"]:
+ if quality_dict[key]["has_good_quality_maxDeviation"]:
+ text_des.append(
+ "The bandOverlaps.lobster file is generated during the LOBSTER run. This "
+ "indicates that the projected wave function is not completely orthonormalized; "
+ "however, the maximal deviation values observed compared to the identity matrix "
+ "is below the threshold of 0.1."
+ )
+ else:
+ text_des.append(
+ "The bandOverlaps.lobster file is generated during the LOBSTER run. This "
+ "indicates that the projected wave function is not completely orthonormalized. "
+ "The maximal deviation value from the identity matrix is {}, and there are "
+ "{} percent k-points above the deviation threshold of 0.1. Please check the "
+ "results of other quality checks like dos comparisons, charges, "
+ "charge spillings before using the results for further "
+ "analysis.".format(
+ quality_dict[key]["max_deviation"],
+ quality_dict[key]["percent_kpoints_abv_limit"],
+ )
+ )
+ else:
+ text_des.append(
+ "The projected wave function is completely orthonormalized as no "
+ "bandOverlaps.lobster file is generated during the LOBSTER run."
+ )
+
+ elif key == "charges":
+ if val:
+ for charge in ["Mulliken", "Loewdin"]:
+ if val[f"BVA_{charge}_agree"]:
+ text_des.append(
+ f"The atomic charge signs from {charge} population analysis agree "
+ "with the bond valence analysis."
+ )
+ if not val[f"BVA_{charge}_agree"]:
+ text_des.append(
+ f"The atomic charge signs from {charge} population analysis do not agree with "
+ "the bond valence analysis."
+ )
+ else:
+ text_des.append(
+ "Oxidation states from BVA analyzer cannot be determined. "
+ "Thus BVA charge comparison is not conducted."
+ )
+
+ elif key == "dos_comparisons":
+ comp_types = []
+ tani_index = []
+ for orb in val:
+ if orb.split("_")[-1] in ["s", "p", "d", "f", "summed"]:
+ comp_types.append(orb.split("_")[-1])
+ tani_index.append(str(val[orb]))
+ text_des.append(
+ "The Tanimoto index from DOS comparisons in the energy range between {}, {} eV "
+ "for {} orbitals are: {}.".format(
+ val["e_range"][0],
+ val["e_range"][1],
+ ", ".join(comp_types),
+ ", ".join(tani_index),
+ )
+ )
+
+ return text_des
+
+[docs] @staticmethod
+ def write_calc_quality_description(calc_quality_text):
+ """Print the calculation quality description to the screen."""
+ print(" ".join(calc_quality_text))
+
+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""This module defines wrapper classes to quickly obtain similarity matrix of input fingerprint objects."""
+
+from __future__ import annotations
+
+import multiprocessing as mp
+import os
+import warnings
+from pathlib import Path
+from typing import NamedTuple
+
+import numpy as np
+import pandas as pd
+from tqdm.autonotebook import tqdm
+
+from lobsterpy.featurize.core import (
+ FeaturizeCharges,
+ FeaturizeCOXX,
+ FeaturizeLobsterpy,
+)
+
+warnings.filterwarnings("ignore")
+
+
+[docs]class BatchSummaryFeaturizer:
+ """
+ Batch Featurizer sets that generates summary features from lobster data.
+
+ Args:
+ path_to_lobster_calcs: path to root directory consisting of all lobster calc
+ path_to_jsons: path to root directory consisting of all lobster lightweight jsons
+ feature_type: set the feature type for moment features.
+ Possible options are "bonding", "antibonding" or "overall"
+ charge_type : set charge type used for computing ionicity. Possible options are
+ "mulliken", "loewdin or "both"
+ bonds: "all_bonds" or "cation_anion_bonds"
+ include_cobi_data : bool stating to include COBICAR.lobster features
+ include_coop_data: bool stating to include COOPCAR.lobster features
+ e_range : range of energy relative to fermi for which moment features needs to be computed
+ n_jobs : parallel processes to run
+
+ Attributes:
+ get_df: A pandas dataframe with summary features
+ """
+
+ def __init__(
+ self,
+ path_to_lobster_calcs: str | Path,
+ path_to_jsons: str | Path | None = None,
+ feature_type: str = "antibonding",
+ charge_type: str = "both",
+ bonds: str = "all",
+ include_cobi_data: bool = False,
+ include_coop_data: bool = False,
+ e_range: list[float] = [-5.0, 0.0],
+ n_jobs: int = 4,
+ ):
+ """
+ Featurize lobster data via multiprocessing for large number of compounds.
+
+ Args:
+ path_to_lobster_calcs: path to root directory consisting of all lobster calc
+ path_to_jsons: path to root directory consisting of all lobster lightweight jsons
+ feature_type: set the feature type for moment features.
+ Possible options are "bonding", "antibonding" or "overall"
+ charge_type : set charge type used for computing ionicity. Possible options are
+ "mulliken", "loewdin or "both"
+ bonds: "all_bonds" or "cation_anion_bonds"
+ include_cobi_data : bool stating to include COBICAR.lobster features
+ include_coop_data: bool stating to include COOPCAR.lobster features
+ e_range : range of energy relative to fermi for which moment features needs to be computed
+ n_jobs : parallel processes to run
+
+ """
+ self.path_to_lobster_calcs = path_to_lobster_calcs
+ self.path_to_jsons = path_to_jsons
+ self.feature_type = feature_type
+ self.charge_type = charge_type
+ self.bonds = bonds
+ self.include_cobi_data = include_cobi_data
+ self.include_coop_data = include_coop_data
+ self.e_range = e_range
+ self.n_jobs = n_jobs
+
+[docs] def _featurizelobsterpy(self, file_name_or_path) -> pd.DataFrame:
+ """
+ Featurize Lobsterpy condensed bonding analysis data.
+
+ if lightweight json file exists loads that or invokes LobsterPy Analysis class
+
+ Returns:
+ A pandas dataframe with ICOHP stats like mean, min, max of relevant bonds and
+ madelung energies
+
+ """
+ if Path(file_name_or_path).is_file():
+ featurize_lobsterpy = FeaturizeLobsterpy(
+ path_to_json=file_name_or_path,
+ bonds=self.bonds,
+ )
+
+ else:
+ featurize_lobsterpy = FeaturizeLobsterpy(
+ path_to_lobster_calc=file_name_or_path,
+ bonds=self.bonds,
+ )
+
+ return featurize_lobsterpy.get_df()
+
+ # return df
+
+[docs] def _featurizecoxx(self, path_to_lobster_calc) -> pd.DataFrame:
+ """
+ Featurize COHP/COBI/COOPCAR data using FeaturizeCOXX.
+
+ Returns:
+ A pandas dataframe with COHP summary stats data mainly weighted ICOHP/ICOOP/ICOBI,
+ Effective interaction number and moment features (center, width, skewness and kurtosis)
+
+ """
+ dir_name = Path(path_to_lobster_calc)
+
+ req_files = {
+ "structure_path": "POSCAR",
+ "coxxcar_path": "COHPCAR.lobster",
+ "icoxxlist_path": "ICOHPLIST.lobster",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ coxxcar_path = req_files.get("coxxcar_path")
+ structure_path = req_files.get("structure_path")
+ icoxxlist_path = req_files.get("icoxxlist_path")
+
+ if (
+ coxxcar_path.exists() # type: ignore
+ and structure_path.exists() # type: ignore
+ and icoxxlist_path.exists() # type: ignore
+ ):
+ coxx = FeaturizeCOXX(
+ path_to_coxxcar=str(coxxcar_path),
+ path_to_icoxxlist=str(icoxxlist_path),
+ path_to_structure=str(structure_path),
+ feature_type=self.feature_type,
+ e_range=self.e_range,
+ )
+
+ df_cohp = coxx.get_summarized_coxx_df()
+ else:
+ raise Exception(
+ "COHPCAR.lobster or POSCAR or ICOHPLIST.lobster file "
+ f"not found in {dir_name.name}"
+ )
+
+ if self.include_cobi_data:
+ req_files = {
+ "coxxcar_path": "COBICAR.lobster",
+ "icoxxlist_path": "ICOBILIST.lobster",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ coxxcar_path = req_files.get("coxxcar_path")
+ icoxxlist_path = req_files.get("icoxxlist_path")
+
+ if coxxcar_path.exists() and icoxxlist_path.exists(): # type: ignore
+ coxx = FeaturizeCOXX(
+ path_to_coxxcar=str(coxxcar_path),
+ path_to_icoxxlist=str(icoxxlist_path),
+ path_to_structure=str(structure_path),
+ feature_type=self.feature_type,
+ e_range=self.e_range,
+ are_cobis=True,
+ )
+
+ df_cobi = coxx.get_summarized_coxx_df()
+
+ else:
+ raise Exception(
+ "COBICAR.lobster or ICOBILIST.lobster file "
+ f"not found in {dir_name.name}"
+ )
+
+ if self.include_coop_data:
+ req_files = {
+ "coxxcar_path": "COOPCAR.lobster",
+ "icoxxlist_path": "ICOOPLIST.lobster",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ coxxcar_path = req_files.get("coxxcar_path")
+ icoxxlist_path = req_files.get("icoxxlist_path")
+
+ if coxxcar_path.exists() and icoxxlist_path.exists(): # type: ignore
+ coxx = FeaturizeCOXX(
+ path_to_coxxcar=str(coxxcar_path),
+ path_to_icoxxlist=str(icoxxlist_path),
+ path_to_structure=str(structure_path),
+ feature_type=self.feature_type,
+ e_range=self.e_range,
+ are_coops=True,
+ )
+
+ df_coop = coxx.get_summarized_coxx_df()
+
+ else:
+ raise Exception(
+ "COOPCAR.lobster or ICOOPLIST.lobster file "
+ f"not found in {dir_name.name}"
+ )
+
+ if self.include_cobi_data and self.include_coop_data:
+ df = pd.concat([df_cohp, df_cobi, df_coop], axis=1)
+ elif self.include_cobi_data and not self.include_coop_data:
+ df = pd.concat([df_cohp, df_cobi], axis=1)
+ elif not self.include_cobi_data and self.include_coop_data:
+ df = pd.concat([df_cohp, df_coop], axis=1)
+ else:
+ df = df_cohp
+
+ return df
+
+[docs] def _featurizecharges(self, path_to_lobster_calc) -> pd.DataFrame:
+ """
+ Featurize CHARGE.lobster.gz data that using FeaturizeCharges.
+
+ Returns:
+ A pandas dataframe with computed ionicity for the structure
+
+ """
+ dir_name = Path(path_to_lobster_calc)
+
+ req_files = {
+ "charge_path": "CHARGE.lobster",
+ "structure_path": "POSCAR",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ charge_path = req_files.get("charge_path")
+ structure_path = req_files.get("structure_path")
+
+ if charge_path.exists() and structure_path.exists(): # type: ignore
+ if self.charge_type == "mulliken":
+ charge_mull = FeaturizeCharges(
+ path_to_charge=str(charge_path),
+ path_to_structure=str(structure_path),
+ charge_type="mulliken",
+ )
+ df = charge_mull.get_df()
+ elif self.charge_type == "loewdin":
+ charge_loew = FeaturizeCharges(
+ path_to_charge=str(charge_path),
+ path_to_structure=str(structure_path),
+ charge_type="loewdin",
+ )
+ df = charge_loew.get_df()
+ elif self.charge_type == "both":
+ charge_mull = FeaturizeCharges(
+ path_to_charge=str(charge_path),
+ path_to_structure=str(structure_path),
+ charge_type="mulliken",
+ )
+ df_mull = charge_mull.get_df()
+
+ charge_loew = FeaturizeCharges(
+ path_to_charge=str(charge_path),
+ path_to_structure=str(structure_path),
+ charge_type="loewdin",
+ )
+ df_loew = charge_loew.get_df()
+
+ df = pd.concat([df_mull, df_loew], axis=1)
+ else:
+ raise Exception(f"CHARGE.lobster or POSCAR not found in {dir_name.name}")
+
+ return df
+
+[docs] def get_df(self) -> pd.DataFrame:
+ """
+ Generate a pandas dataframe with summary features extracted from LOBSTER files.
+
+ Uses multiprocessing to speed up the process.
+
+ Returns:
+ Returns a pandas dataframe
+
+ """
+ if self.path_to_jsons:
+ file_name_or_path = [
+ os.path.join(self.path_to_jsons, f)
+ for f in os.listdir(self.path_to_jsons)
+ if not f.startswith("t")
+ and not f.startswith(".")
+ and not os.path.isdir(f)
+ ]
+
+ elif self.path_to_lobster_calcs and not self.path_to_jsons:
+ file_name_or_path = [
+ os.path.join(self.path_to_lobster_calcs, f)
+ for f in os.listdir(self.path_to_lobster_calcs)
+ if not f.startswith("t")
+ and not f.startswith(".")
+ and os.path.isdir(os.path.join(self.path_to_lobster_calcs, f))
+ ]
+
+ with mp.Pool(processes=self.n_jobs, maxtasksperchild=1) as pool:
+ results = tqdm(
+ pool.imap_unordered(
+ self._featurizelobsterpy, file_name_or_path, chunksize=1
+ ),
+ total=len(file_name_or_path),
+ desc="Generating LobsterPy summary stats",
+ )
+ pool.close()
+ pool.join()
+ row = []
+ for result in results:
+ row.append(result) # noqa: PERF402
+
+ df_lobsterpy = pd.concat(row)
+ df_lobsterpy.sort_index(inplace=True) # noqa: PD002
+
+ paths = [
+ os.path.join(self.path_to_lobster_calcs, f)
+ for f in os.listdir(self.path_to_lobster_calcs)
+ if not f.startswith("t")
+ and not f.startswith(".")
+ and os.path.isdir(os.path.join(self.path_to_lobster_calcs, f))
+ ]
+
+ with mp.Pool(processes=self.n_jobs, maxtasksperchild=1) as pool:
+ results = tqdm(
+ pool.imap_unordered(self._featurizecoxx, paths, chunksize=1),
+ total=len(paths),
+ desc="Generating COHP/COOP/COBI summary stats",
+ )
+ pool.close()
+ pool.join()
+ row = []
+ for result in results:
+ row.append(result)
+
+ df_coxx = pd.concat(row)
+ df_coxx.sort_index(inplace=True) # noqa: PD002
+
+ with mp.Pool(processes=self.n_jobs, maxtasksperchild=1) as pool:
+ results = tqdm(
+ pool.imap_unordered(self._featurizecharges, paths, chunksize=1),
+ total=len(paths),
+ desc="Generating charge based features",
+ )
+ pool.close()
+ pool.join()
+ row = []
+ for result in results:
+ row.append(result)
+
+ df_charges = pd.concat(row)
+ df_charges.sort_index(inplace=True) # noqa: PD002
+
+ return pd.concat([df_lobsterpy, df_coxx, df_charges], axis=1)
+
+ # return df
+
+
+[docs]class BatchCoxxFingerprint:
+ """
+ BatchFeaturizer to generate COHP/COOP/COBI fingerprints and Tanimoto index similarity matrix.
+
+ Args:
+ path_to_lobster_calcs: path to root directory consisting of all lobster calc
+ feature_type: set the feature type for moment features.
+ Possible options are "bonding", "antibonding" or "overall"
+ label_list: bond labels list for which fingerprints needs to be generated.
+ tanimoto : bool to state to compute tanimoto index between fingerprint objects
+ normalize: bool to state to normalize the fingerprint data
+ n_bins: sets number for bins for fingerprint objects
+ e_range : range of energy relative to fermi for which moment features needs to be computed
+ n_jobs : number of parallel processes to run
+ fingerprint_for: Possible options are 'cohp/cobi/coop'.
+ Based on this fingerprints will be computed for COHPCAR/COOBICAR/COOPCAR.lobster files
+
+ Attributes:
+ fingerprint_df: A pandas dataframe with fingerprint objects
+ get_similarity_matrix_df: A symmetric pandas dataframe consisting of
+ similarity index (tanimoto/normalized dot product/dot product)
+ computed between all pairs of compounds
+ """
+
+ def __init__(
+ self,
+ path_to_lobster_calcs: str | Path,
+ feature_type: str = "overall",
+ label_list: list[str] | None = None,
+ tanimoto: bool = True,
+ normalize: bool = True,
+ spin_type: str = "summed",
+ n_bins: int = 56,
+ e_range: list[float] = [-15.0, 0.0],
+ n_jobs=4,
+ fingerprint_for: str = "cohp",
+ ):
+ """
+ Generate COHP/COOP/COBI fingerprints and pair-wise Tanimoto index similarity matrix.
+
+ Args:
+ path_to_lobster_calcs: path to root directory consisting of all lobster calc
+ feature_type: set the feature type for moment features.
+ Possible options are "bonding", "antibonding" or "overall"
+ label_list: bond labels list for which fingerprints needs to be generated.
+ tanimoto : bool to state to compute tanimoto index between fingerprint objects
+ normalize: bool to state to normalize the fingerprint data
+ spin_type: can be summed or up or down.
+ n_bins: sets number for bins for fingerprint objects
+ e_range : range of energy relative to fermi for which moment features needs to be computed
+ n_jobs : number of parallel processes to run
+ fingerprint_for: Possible options are 'cohp/cobi/coop'.
+ Based on this fingerprints will be computed for COHPCAR/COOBICAR/COOPCAR.lobster files
+
+ """
+ self.path_to_lobster_calcs = path_to_lobster_calcs
+ self.feature_type = feature_type
+ self.tanimoto = tanimoto
+ self.normalize = normalize
+ self.label_list = label_list
+ self.spin_type = spin_type
+ self.n_bins = n_bins
+ self.e_range = e_range
+ self.n_jobs = n_jobs
+ self.fingerprint_for = fingerprint_for
+
+ self.fingerprint_df = self._get_fingerprints_df()
+
+[docs] def get_similarity_matrix_df(self) -> pd.DataFrame:
+ """
+ Compute pairwise similarity index for each fingerprint object in input dataframe.
+
+ Returns:
+ A Pandas dataframe
+ """
+ matrix = np.full(
+ (self.fingerprint_df.shape[0], self.fingerprint_df.shape[0]), np.nan
+ )
+ for i, (_, col) in enumerate(self.fingerprint_df.iterrows()):
+ for j, (_, col1) in enumerate(self.fingerprint_df.iterrows()):
+ if self.tanimoto:
+ simi = self._get_fp_similarity(
+ col["COXX_FP"],
+ col1["COXX_FP"],
+ tanimoto=self.tanimoto,
+ normalize=False,
+ )
+ else:
+ simi = self._get_fp_similarity(
+ col["COXX_FP"],
+ col1["COXX_FP"],
+ tanimoto=False,
+ normalize=True,
+ )
+ matrix[i][j] = simi
+
+ return pd.DataFrame(
+ matrix,
+ index=list(self.fingerprint_df.index),
+ columns=list(self.fingerprint_df.index),
+ )
+
+[docs] @staticmethod
+ def _fp_to_dict(fp) -> dict:
+ """
+ Convert a fingerprint obj into a dictionary.
+
+ Args:
+ fp: The fingerprint to be converted into a dictionary
+
+ Returns:
+ dict: A dict of the fingerprint Keys=type, Values=np.ndarray(energies, cohp)
+ """
+ fp_dict = {}
+ fp_dict[fp[2]] = np.array([fp[0], fp[1]], dtype="object").T
+
+ return fp_dict
+
+[docs] @staticmethod
+ def _get_fp_similarity(
+ fp1: NamedTuple,
+ fp2: NamedTuple,
+ col: int = 1,
+ pt: int | str = "All",
+ normalize: bool = False,
+ tanimoto: bool = True,
+ ) -> float:
+ """
+ Calculate the similarity index (dot product) of two fingerprints.
+
+ Args:
+ fp1 (NamedTuple): The 1st dos fingerprint object
+ fp2 (NamedTuple): The 2nd dos fingerprint object
+ col (int): The item in the fingerprints (0:energies,1: coxxs) to take the dot product of (default is 1)
+ pt (int or str) : The index of the point that the dot product is to be taken (default is All)
+ normalize (bool): If True normalize the scalar product to 1 (default is False)
+ tanimoto (bool): If True will compute Tanimoto index (default is False)
+
+ Raises:
+ ValueError: If both tanimoto and normalize are set to True.
+
+ Returns:
+ Similarity index (float): The value of dot product
+
+ """
+ fp1_dict = (
+ BatchCoxxFingerprint._fp_to_dict(fp1) if not isinstance(fp1, dict) else fp1
+ )
+
+ fp2_dict = (
+ BatchCoxxFingerprint._fp_to_dict(fp2) if not isinstance(fp2, dict) else fp2
+ )
+
+ if pt == "All":
+ vec1 = np.array([pt[col] for pt in fp1_dict.values()]).flatten()
+ vec2 = np.array([pt[col] for pt in fp2_dict.values()]).flatten()
+ else:
+ vec1 = fp1_dict[fp1[2][pt]][col]
+ vec2 = fp2_dict[fp2[2][pt]][col]
+
+ if not normalize and tanimoto:
+ rescale = (
+ np.linalg.norm(vec1) ** 2
+ + np.linalg.norm(vec2) ** 2
+ - np.dot(vec1, vec2)
+ )
+
+ elif not tanimoto and normalize:
+ rescale = np.linalg.norm(vec1) * np.linalg.norm(vec2)
+
+ elif not tanimoto and not normalize:
+ rescale = 1.0
+
+ else:
+ raise ValueError(
+ "Cannot compute similarity index. Please set either normalize=True or tanimoto=True or both to False."
+ )
+ return np.dot(vec1, vec2) / rescale
+
+[docs] def _fingerprint_df(self, path_to_lobster_calc) -> pd.DataFrame:
+ """
+ Get fingerprint object dataframe via FeaturizeCOXX.get_coxx_fingerprint_df.
+
+ Also helps to generate the data used for fingerprint generation
+
+ Returns:
+ A pandas dataframe with COXX fingerprint object
+
+ """
+ dir_name = Path(path_to_lobster_calc)
+
+ if self.fingerprint_for.upper() == "COBI":
+ req_files = {
+ "coxxcar_path": "COBICAR.lobster",
+ "icoxxlist_path": "ICOBILIST.lobster",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ coxxcar_path = req_files.get("coxxcar_path")
+ icoxxlist_path = req_files.get("icoxxlist_path")
+ are_cobis = True
+ are_coops = False
+
+ elif self.fingerprint_for.upper() == "COOP":
+ req_files = {
+ "coxxcar_path": "COOPCAR.lobster",
+ "icoxxlist_path": "ICOOPLIST.lobster",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ coxxcar_path = req_files.get("coxxcar_path")
+ icoxxlist_path = req_files.get("icoxxlist_path")
+ are_cobis = False
+ are_coops = True
+
+ else:
+ req_files = {
+ "coxxcar_path": "COHPCAR.lobster",
+ "icoxxlist_path": "ICOHPLIST.lobster",
+ }
+ for file, default_value in req_files.items():
+ file_path = dir_name / default_value
+ req_files[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files[file] = gz_file_path # type: ignore
+
+ coxxcar_path = req_files.get("coxxcar_path")
+ icoxxlist_path = req_files.get("icoxxlist_path")
+ are_cobis = False
+ are_coops = False
+
+ structure_path = dir_name / "POSCAR"
+ if not structure_path.exists():
+ gz_file_path = structure_path.with_name(structure_path.name + ".gz")
+ if gz_file_path.exists():
+ structure_path = gz_file_path
+
+ coxx = FeaturizeCOXX(
+ path_to_coxxcar=str(coxxcar_path),
+ path_to_icoxxlist=str(icoxxlist_path),
+ path_to_structure=str(structure_path),
+ feature_type=self.feature_type,
+ e_range=self.e_range,
+ are_coops=are_coops,
+ are_cobis=are_cobis,
+ )
+
+ return coxx.get_coxx_fingerprint_df(
+ spin_type=self.spin_type,
+ n_bins=self.n_bins,
+ normalize=self.normalize,
+ label_list=self.label_list,
+ )
+
+[docs] def _get_fingerprints_df(self) -> pd.DataFrame:
+ """
+ Generate fingerprint objects dataframe using BatchCoxxFingerprint._fingerprint_df.
+
+ Returns:
+ A pandas dataframe with COXX fingerprint objects
+
+ """
+ paths = [
+ os.path.join(self.path_to_lobster_calcs, f)
+ for f in os.listdir(self.path_to_lobster_calcs)
+ if not f.startswith("t")
+ and not f.startswith(".")
+ and os.path.isdir(os.path.join(self.path_to_lobster_calcs, f))
+ ]
+
+ with mp.Pool(processes=self.n_jobs, maxtasksperchild=1) as pool:
+ results = tqdm(
+ pool.imap_unordered(self._fingerprint_df, paths, chunksize=1),
+ total=len(paths),
+ desc=f"Generating {self.fingerprint_for.upper()} fingerprints",
+ )
+ pool.close()
+ pool.join()
+ row = []
+ for result in results:
+ row.append(result) # noqa: PERF402
+
+ df = pd.concat(row)
+ df.sort_index(inplace=True) # noqa: PD002
+
+ return df
+
+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""This module defines classes to featurize Lobster data ready to be used for ML studies."""
+
+from __future__ import annotations
+
+import gzip
+import json
+import warnings
+from collections import namedtuple
+from pathlib import Path
+
+import numpy as np
+import numpy.typing as npt
+import pandas as pd
+from mendeleev import element
+from pymatgen.core.structure import Structure
+from pymatgen.electronic_structure.cohp import CompleteCohp
+from pymatgen.electronic_structure.core import Spin
+from pymatgen.io.lobster import Charge, Icohplist, MadelungEnergies
+from scipy.integrate import trapezoid
+
+from lobsterpy.cohp.analyze import Analysis
+
+warnings.filterwarnings("ignore")
+
+
+[docs]class FeaturizeLobsterpy:
+ """
+ Class to featurize lobsterpy data.
+
+ Args:
+ path_to_lobster_calc: path to parent directory containing lobster calc outputs
+ path_to_json: path to lobster lightweight json
+ bonds: "all" or "cation-anion" bonds
+ Attributes:
+ get_df: returns a pandas dataframe with relevant icohp statistical data as columns from
+ lobsterpy automatic bonding analysis
+
+ """
+
+ def __init__(
+ self,
+ path_to_lobster_calc: str | Path | None = None,
+ path_to_json: str | Path | None = None,
+ bonds: str = "all",
+ ):
+ """
+ Extract features from Lobster calculations.
+
+ Args:
+ path_to_lobster_calc: path to parent directory containing lobster calc outputs
+ path_to_json: path to lobster lightweight json
+ bonds: "all" or "cation-anion" bonds
+
+ """
+ self.path_to_json = path_to_json
+ self.path_to_lobster_calc = path_to_lobster_calc
+ self.bonds = bonds
+
+[docs] def get_df(self, ids: str | None = None) -> pd.DataFrame:
+ """
+ Featurize LobsterPy condensed bonding analysis data.
+
+ Returns:
+ Returns a pandas dataframe with lobsterpy icohp statistics
+
+ """
+ if self.path_to_json and not self.path_to_lobster_calc:
+ # read the lightweight lobster json files using read_lobster_lightweight_json method
+ data = FeaturizeLobsterpy.read_lobster_lightweight_json(
+ path_to_json=self.path_to_json
+ )
+ if not ids:
+ ids = Path(self.path_to_json).name.split(".")[0]
+
+ elif self.path_to_lobster_calc and not self.path_to_json:
+ # get lobsterpy condensed bonding analysis data using get_lobsterpy_cba_dict method
+ data = FeaturizeLobsterpy.get_lobsterpy_cba_dict(
+ path_to_lobster_calc=self.path_to_lobster_calc, bonds=self.bonds
+ )
+
+ if not ids:
+ ids = Path(self.path_to_lobster_calc).name
+
+ else:
+ raise ValueError(
+ "Please provide either path to lightweight lobster jsons or path to lobster calc"
+ )
+ # define a pandas dataframe
+ df = pd.DataFrame(index=[ids])
+
+ icohp_mean = []
+ icohp_sum = []
+ bond = []
+ antibond = []
+
+ # extract lobsterpy icohp related data for bond type specified
+ # Results will differ for "all" and "cation-anion" mode.
+ # In "all" bonds mode, the bonds will come up twice, also
+ # cation-cation, anion-anion bonds will also be considered
+
+ if self.bonds == "all":
+ bond_type = "all_bonds"
+ elif self.bonds == "cation-anion":
+ bond_type = "cation_anion_bonds"
+
+ if (
+ not data[bond_type]["lobsterpy_data"]
+ or not data[bond_type]["lobsterpy_data"]["sites"]
+ ):
+ raise Exception(
+ f"No {self.bonds} bonds detected for {ids} structure. "
+ "Please switch to `all` bonds mode"
+ )
+
+ for site_data in data[bond_type]["lobsterpy_data"]["sites"].values():
+ if site_data["bonds"]:
+ for bond_data in site_data["bonds"].values():
+ icohp_mean.append(float(bond_data["ICOHP_mean"]))
+ icohp_sum.append(float(bond_data["ICOHP_sum"]))
+ bond.append(bond_data["bonding"]["perc"])
+ antibond.append(bond_data["antibonding"]["perc"])
+
+ # add ICOHP stats data (mean, min, max, standard deviation) as columns to the dataframe
+ df.loc[ids, "Icohp_mean_avg"] = np.mean(icohp_mean)
+ df.loc[ids, "Icohp_mean_max"] = np.max(icohp_mean)
+ df.loc[ids, "Icohp_mean_min"] = np.min(icohp_mean)
+ df.loc[ids, "Icohp_mean_std"] = np.std(icohp_mean)
+
+ df.loc[ids, "Icohp_sum_avg"] = np.mean(icohp_sum)
+ df.loc[ids, "Icohp_sum_max"] = np.max(icohp_sum)
+ df.loc[ids, "Icohp_sum_min"] = np.min(icohp_sum)
+ df.loc[ids, "Icohp_sum_std"] = np.std(icohp_sum)
+
+ df.loc[ids, "bonding_perc_avg"] = np.mean(bond)
+ df.loc[ids, "bonding_perc_max"] = np.max(bond)
+ df.loc[ids, "bonding_perc_min"] = np.min(bond)
+ df.loc[ids, "bonding_perc_std"] = np.std(bond)
+
+ df.loc[ids, "antibonding_perc_avg"] = np.mean(antibond)
+ df.loc[ids, "antibonding_perc_min"] = np.min(antibond)
+ df.loc[ids, "antibonding_perc_max"] = np.max(antibond)
+ df.loc[ids, "antibonding_perc_std"] = np.std(antibond)
+
+ # add madelung energies for the structure
+ df.loc[ids, "Madelung_Mull"] = data["madelung_energies"]["Mulliken"]
+ df.loc[ids, "Madelung_Loew"] = data["madelung_energies"]["Loewdin"]
+
+ return df
+
+[docs] @staticmethod
+ def read_lobster_lightweight_json(path_to_json: str | Path) -> dict:
+ """
+ Read the lightweight JSON.gz files and return a Python dictionary object.
+
+ Args:
+ path_to_json: path to lobsterpy lightweight json file
+
+ Returns:
+ Returns a dictionary with lobster summarized bonding analysis data
+
+ """
+ with gzip.open(str(path_to_json), "rb") as f:
+ data = json.loads(f.read().decode("utf-8"))
+
+ lobster_data = {}
+ for item in data:
+ lobster_data.update(item)
+
+ return lobster_data
+
+[docs] @staticmethod
+ def get_lobsterpy_cba_dict(path_to_lobster_calc: str | Path, bonds: str) -> dict:
+ """
+ Generate a Python dictionary object using the Analysis class with condensed bonding analysis data.
+
+ Args:
+ path_to_lobster_calc: path to lobsterpy lightweight json file
+ bonds: "all" or "cation-anion" bonds
+
+ Returns:
+ Returns a dictionary with lobster summarized bonding analysis data
+
+ """
+ dir_name = Path(str(path_to_lobster_calc))
+
+ # check if files are compressed (.gz) and update file paths
+ req_files_lobsterpy = {
+ "structure_path": "POSCAR",
+ "cohpcar_path": "COHPCAR.lobster",
+ "icohplist_path": "ICOHPLIST.lobster",
+ "charge_path": "CHARGE.lobster",
+ }
+
+ for file, default_value in req_files_lobsterpy.items():
+ file_path = dir_name / default_value
+ req_files_lobsterpy[file] = file_path # type: ignore
+ if not file_path.exists():
+ gz_file_path = file_path.with_name(file_path.name + ".gz")
+ if gz_file_path.exists():
+ req_files_lobsterpy[file] = gz_file_path # type: ignore
+ else:
+ raise Exception(
+ "Path provided for Lobster calc directory seems incorrect."
+ "It does not contain COHPCAR.lobster, ICOHPLIST.lobster, POSCAR and "
+ "CHARGE.lobster files needed for automatic analysis using LobsterPy"
+ )
+
+ cohpcar_path = req_files_lobsterpy.get("cohpcar_path")
+ charge_path = req_files_lobsterpy.get("charge_path")
+ structure_path = req_files_lobsterpy.get("structure_path")
+ icohplist_path = req_files_lobsterpy.get("icohplist_path")
+
+ which_bonds = bonds.replace("-", "_")
+ bond_type = f"{which_bonds}_bonds"
+
+ try:
+ analyse = Analysis(
+ path_to_poscar=str(structure_path),
+ path_to_icohplist=str(icohplist_path),
+ path_to_cohpcar=str(cohpcar_path),
+ path_to_charge=str(charge_path),
+ summed_spins=False, # we will always use spin polarization here
+ cutoff_icohp=0.10,
+ which_bonds=which_bonds,
+ )
+
+ data = {bond_type: {"lobsterpy_data": analyse.condensed_bonding_analysis}}
+ except ValueError:
+ data = {bond_type: {"lobsterpy_data": {}}}
+
+ madelung_energies_path = dir_name / "MadelungEnergies.lobster"
+ # check if .gz file exists and update Madelung Energies path
+ if not madelung_energies_path.exists():
+ gz_file_path = madelung_energies_path.with_name(
+ madelung_energies_path.name + ".gz"
+ )
+ if gz_file_path.exists():
+ madelung_energies_path = gz_file_path
+
+ if madelung_energies_path.exists():
+ madelung_obj = MadelungEnergies(filename=str(madelung_energies_path))
+
+ madelung_energies = {
+ "Mulliken": madelung_obj.madelungenergies_Mulliken,
+ "Loewdin": madelung_obj.madelungenergies_Loewdin,
+ "Ewald_splitting": madelung_obj.ewald_splitting,
+ }
+ data["madelung_energies"] = madelung_energies
+
+ else:
+ warnings.warn(
+ "MadelungEnergies.lobster file not found in Lobster calc directory provided"
+ " Will set Madelung Engeries for crystal structure values to NaN"
+ )
+ madelung_energies = {
+ "Mulliken": np.nan,
+ "Loewdin": np.nan,
+ "Ewald_splitting": np.nan,
+ }
+
+ data["madelung_energies"] = madelung_energies
+
+ return data
+
+
+coxx_fingerprint = namedtuple(
+ "coxx_fingerprint", "energies coxx fp_type spin_type n_bins bin_width"
+)
+
+
+[docs]class FeaturizeCOXX:
+ """
+ Class to featurize COHPCAR, COBICAR or COOPCAR data.
+
+ Args:
+ path_to_coxxcar: path to COXXCAR.lobster (e.g., "COXXCAR.lobster")
+ path_to_icoxxlist : path to ICOXXLIST.lobster (e.g., "ICOXXLIST.lobster")
+ path_to_structure : path to structure file (e.g., "POSCAR")
+ feature_type: set the feature type for moment features and fingerprints.
+ Possible options are "bonding", "antibonding" or "overall"
+ are_cobis : bool indicating if file contains COBI/ICOBI data
+ are_coops : bool indicating if file contains COOP/ICOOP data
+ e_range : range of energy relative to fermi for which moment features needs to be computed
+
+ Attributes:
+ get_df: pandas dataframe
+ get_coxx_fingerprint_df: pandas dataframe
+
+ """
+
+ def __init__(
+ self,
+ path_to_coxxcar: str | Path,
+ path_to_icoxxlist: str | Path,
+ path_to_structure: str | Path,
+ feature_type: str,
+ e_range: list[float] = [-10.0, 0.0],
+ are_cobis: bool = False,
+ are_coops: bool = False,
+ ):
+ """
+ Featurize COHPCAR, COBICAR or COOPCAR data.
+
+ Args:
+ path_to_coxxcar: path to COXXCAR.lobster (e.g., "COXXCAR.lobster")
+ path_to_icoxxlist : path to ICOXXLIST.lobster (e.g., "ICOXXLIST.lobster")
+ path_to_structure : path to structure file (e.g., "POSCAR")
+ feature_type: set the feature type for moment features and fingerprints.
+ Possible options are "bonding", "antibonding" or "overall"
+ are_cobis : bool indicating if file contains COBI/ICOBI data
+ are_coops : bool indicating if file contains COOP/ICOOP data
+ e_range : range of energy relative to fermi for which moment features needs to be computed
+
+ """
+ self.path_to_coxxcar = path_to_coxxcar
+ self.path_to_icoxxlist = path_to_icoxxlist
+ self.path_to_structure = path_to_structure
+ self.feature_type = feature_type
+ self.e_range = e_range
+ self.are_cobis = are_cobis
+ self.are_coops = are_coops
+ self.icoxxlist = Icohplist(
+ filename=self.path_to_icoxxlist,
+ are_cobis=self.are_cobis,
+ are_coops=self.are_coops,
+ )
+ self.completecoxx = CompleteCohp.from_file(
+ fmt="LOBSTER",
+ filename=self.path_to_coxxcar,
+ structure_file=self.path_to_structure,
+ are_cobis=self.are_cobis,
+ are_coops=self.are_coops,
+ )
+
+[docs] def get_coxx_fingerprint_df(
+ self,
+ ids: str | None = None,
+ label_list: list[str] | None = None,
+ per_bond: bool = True,
+ spin_type: str = "summed",
+ binning: bool = True,
+ n_bins: int = 56,
+ normalize: bool = True,
+ ) -> pd.DataFrame:
+ """
+ Generate a COXX fingerprints dataframe.
+
+ Args:
+ ids: sets index of pandas dataframe
+ spin_type: Specify spin type. Can accept '{summed/up/down}'
+ (default is summed)
+ binning: If true coxxs will be binned
+ n_bins: Number of bins to be used in the fingerprint (default is 256)
+ normalize: If true, normalizes the area under fp to equal to 1 (default is True)
+ label_list: Specify bond labels as a list for which cohp fingerprints are needed
+ per_bond: Will scale cohp values by number of bonds i.e. length of label_list arg
+ (Only affects when label_list is not None)
+
+
+ Raises:
+ Exception: If spin_type is not one of the accepted values {summed/up/down}.
+ ValueError: If feature_type is not one of the accepted values {bonding/antibonding/overall}.
+
+ Returns:
+ A pandas dataframe with the COXX fingerprint
+ of format (energies, coxx, fp_type, spin_type, n_bins, bin_width)
+
+ """
+ coxxcar_obj = self.completecoxx
+
+ energies = coxxcar_obj.energies - coxxcar_obj.efermi
+
+ min_e = self.e_range[0]
+ max_e = self.e_range[-1]
+
+ if max_e is None:
+ max_e = np.max(energies)
+
+ if min_e is None:
+ min_e = np.min(energies)
+
+ if label_list:
+ divisor = len(label_list) if per_bond else 1
+ coxxcar_obj = coxxcar_obj.get_summed_cohp_by_label_list(
+ label_list, divisor=divisor
+ ).get_cohp()
+ else:
+ coxxcar_obj = coxxcar_obj.get_cohp()
+
+ if spin_type == "up":
+ coxx_all = coxxcar_obj[Spin.up]
+ elif spin_type == "down":
+ if Spin.down in coxxcar_obj:
+ coxx_all = coxxcar_obj[Spin.down]
+ else:
+ raise ValueError(
+ "LOBSTER calculation is non-spin polarized. Please switch spin_type to `up`"
+ )
+ elif spin_type == "summed":
+ if Spin.down in coxxcar_obj:
+ coxx_all = coxxcar_obj[Spin.up] + coxxcar_obj[Spin.down]
+ else:
+ coxx_all = coxxcar_obj[Spin.up]
+ else:
+ raise Exception(
+ "Check the spin_type argument.Possible options are summed/up/down"
+ )
+
+ coxx_dict = {}
+ tol = 1e-6
+ if not self.are_cobis and not self.are_coops:
+ coxx_dict["bonding"] = np.array(
+ [scohp if scohp <= tol else 0 for scohp in coxx_all]
+ )
+ coxx_dict["antibonding"] = np.array(
+ [scohp if scohp >= tol else 0 for scohp in coxx_all]
+ )
+ else:
+ coxx_dict["antibonding"] = np.array(
+ [scohp if scohp <= tol else 0 for scohp in coxx_all]
+ )
+ coxx_dict["bonding"] = np.array(
+ [scohp if scohp >= tol else 0 for scohp in coxx_all]
+ )
+
+ coxx_dict["overall"] = coxx_all
+
+ try:
+ if ids:
+ df = pd.DataFrame(index=[ids], columns=["COXX_FP"])
+ else:
+ ids = Path(self.path_to_coxxcar).parent.name
+ df = pd.DataFrame(index=[ids], columns=["COXX_FP"])
+
+ coxxs = coxx_dict[self.feature_type]
+ if len(energies) < n_bins:
+ inds = np.where((energies >= min_e - tol) & (energies <= max_e + tol))
+ fp = coxx_fingerprint(
+ energies[inds],
+ coxxs[inds],
+ self.feature_type,
+ spin_type,
+ len(energies),
+ np.diff(energies)[0],
+ )
+
+ df.at[ids, "COXX_FP"] = fp # noqa: PD008
+
+ return df
+
+ if binning:
+ ener_bounds = np.linspace(min_e, max_e, n_bins + 1)
+ ener = ener_bounds[:-1] + (ener_bounds[1] - ener_bounds[0]) / 2.0
+ bin_width = np.diff(ener)[0]
+ else:
+ ener_bounds = np.array(energies)
+ ener = np.append(energies, [energies[-1] + np.abs(energies[-1]) / 10])
+ n_bins = len(energies)
+ bin_width = np.diff(energies)[0]
+
+ coxx_rebin = np.zeros(ener.shape)
+
+ for ii, e1, e2 in zip(range(len(ener)), ener_bounds[0:-1], ener_bounds[1:]):
+ inds = np.where((energies >= e1) & (energies < e2))
+ coxx_rebin[ii] = np.sum(coxxs[inds])
+ if normalize: # scale DOS bins to make area under histogram equal 1
+ area = np.sum([abs(coxx) * bin_width for coxx in coxx_rebin])
+ coxx_rebin_sc = coxx_rebin / area
+ else:
+ coxx_rebin_sc = coxx_rebin
+
+ fp = coxx_fingerprint(
+ np.array([ener]),
+ coxx_rebin_sc,
+ self.feature_type,
+ spin_type,
+ n_bins,
+ bin_width,
+ )
+
+ df.at[ids, "COXX_FP"] = fp # noqa: PD008
+
+ return df
+
+ except KeyError:
+ raise ValueError(
+ "Please recheck fp_type requested argument.Possible options are bonding/antibonding/overall"
+ )
+
+[docs] def _calculate_wicoxx_ein(self) -> tuple[float, float, float, float]:
+ r"""
+ Calculate weighted ICOXX (xx ==hp,op,bi) and EIN of crystal structure based on work by Nelson et al.
+
+ More details can be found here: https://doi.org/10.1016/B978-0-12-823144-9.00120-5
+
+ \\[\bar{COHP}(E) = \\sum_{i} \\left( w_{i} \\cdot COHP_{i}(E) \right)\\]
+
+ where:
+ - \\(w_{i} = \frac{ICOHP_{i}}{ICOHP_{\text{total}}}\\)
+
+ \\[EIN = \frac{{ICOHP_{\text{{total}}}}}{{\\overline{ICOHP}}} \\cdot \frac{2}{{N_{\text{{Atoms}}}}}\\]
+
+ where:
+ - \\(\\overline{ICOHP}\\) represents the average ICOHP value.
+ - \\(ICOHP_{\text{{total}}}\\) is the total ICOHP value.
+ - \\(N_{\text{{Atoms}}}\\) is the number of atoms in the system.
+
+ Returns:
+ Percent bonding, Percent anti-bonding, weighted icoxx, effective interaction number
+
+ """
+ list_labels = list(self.icoxxlist.icohplist.keys())
+ # Compute sum of icohps
+ icoxx_total = self.icoxxlist.icohpcollection.get_summed_icohp_by_label_list(
+ list_labels
+ )
+
+ summed_weighted_coxx = []
+
+ for lab in list_labels:
+ for cohp in self.completecoxx.get_cohp_by_label(
+ f"{lab}", summed_spin_channels=True
+ ).cohp.values():
+ coxx = cohp
+ icoxx = self.icoxxlist.icohpcollection.get_icohp_by_label(
+ lab, summed_spin_channels=True
+ )
+ weight = (
+ icoxx / icoxx_total
+ ) # calculate the weights based on icohp contri to total icohp of the structure
+ weighted_coxx = weight * coxx
+ summed_weighted_coxx.append(weighted_coxx)
+
+ summed = np.sum(summed_weighted_coxx, axis=0)
+ # below fermi
+ indices = self.completecoxx.energies <= self.completecoxx.efermi
+ en_bf = self.completecoxx.energies[indices]
+ coxx_bf = summed[indices]
+
+ w_icoxx = trapezoid(coxx_bf, en_bf)
+
+ ein = (icoxx_total / w_icoxx) * (
+ 2 / self.completecoxx.structure.num_sites
+ ) # calc effective interaction number
+
+ # percent bonding-anitbonding
+ tol = 1e-6
+ if not self.icoxxlist.are_cobis and not self.icoxxlist.are_coops:
+ bonding_indices = coxx_bf <= tol
+ antibonding_indices = coxx_bf >= tol
+ bnd = abs(trapezoid(en_bf[bonding_indices], coxx_bf[bonding_indices]))
+ antibnd = abs(
+ trapezoid(en_bf[antibonding_indices], coxx_bf[antibonding_indices])
+ )
+ per_bnd = (bnd / (bnd + antibnd)) * 100
+ per_antibnd = (antibnd / (bnd + antibnd)) * 100
+
+ elif self.icoxxlist.are_cobis or self.icoxxlist.are_coops:
+ bonding_indices = coxx_bf >= tol
+ antibonding_indices = coxx_bf <= tol
+ bnd = abs(trapezoid(coxx_bf[bonding_indices], en_bf[bonding_indices]))
+ antibnd = abs(
+ trapezoid(coxx_bf[antibonding_indices], en_bf[antibonding_indices])
+ )
+ per_bnd = (bnd / (bnd + antibnd)) * 100
+ per_antibnd = (antibnd / (bnd + antibnd)) * 100
+
+ return per_bnd, per_antibnd, w_icoxx, ein
+
+[docs] def _calc_moment_features(
+ self, label_list: list[str] | None = None, per_bond=True
+ ) -> tuple[float, float, float, float]:
+ """
+ Calculate band center,width, skewness, and kurtosis of the COXX.
+
+ Args:
+ label_list: List of bond labels
+ per_bond: Will scale cohp values by number of bonds i.e .length of label_list arg
+ (Only affects when label_list is not None)
+
+ Returns:
+ coxx center,width, skewness, and kurtosis in eV
+ """
+ if label_list:
+ divisor = len(label_list) if per_bond else 1
+ coxxcar = self.completecoxx.get_summed_cohp_by_label_list(
+ label_list, divisor=divisor
+ ).get_cohp()
+
+ else:
+ coxxcar = self.completecoxx.get_cohp()
+
+ coxx_all = (
+ coxxcar[Spin.up] + coxxcar[Spin.down]
+ if Spin.down in coxxcar
+ else coxxcar[Spin.up]
+ )
+ energies = self.completecoxx.energies - self.completecoxx.efermi
+
+ coxx_dict = {}
+ tol = 1e-6
+ if not self.are_cobis and not self.are_coops:
+ coxx_dict["bonding"] = np.array(
+ [scohp if scohp <= tol else 0 for scohp in coxx_all]
+ )
+ coxx_dict["antibonding"] = np.array(
+ [scohp if scohp >= tol else 0 for scohp in coxx_all]
+ )
+ coxx_dict["overall"] = coxx_all
+ else:
+ coxx_dict["antibonding"] = np.array(
+ [scohp if scohp <= tol else 0 for scohp in coxx_all]
+ )
+ coxx_dict["bonding"] = np.array(
+ [scohp if scohp >= tol else 0 for scohp in coxx_all]
+ )
+ coxx_dict["overall"] = coxx_all
+
+ try:
+ coxx_center = self._get_coxx_center(
+ coxx=coxx_dict[self.feature_type],
+ energies=energies,
+ e_range=self.e_range,
+ )
+ coxx_width = np.sqrt(
+ self._get_n_moment(
+ n=2,
+ coxx=coxx_dict[self.feature_type],
+ energies=energies,
+ e_range=self.e_range,
+ )
+ )
+ coxx_skew = self._get_n_moment(
+ n=3,
+ coxx=coxx_dict[self.feature_type],
+ energies=energies,
+ e_range=self.e_range,
+ ) / self._get_n_moment(
+ n=2,
+ coxx=coxx_dict[self.feature_type],
+ energies=energies,
+ e_range=self.e_range,
+ ) ** (
+ 3 / 2
+ )
+
+ coxx_kurt = (
+ self._get_n_moment(
+ n=4,
+ coxx=coxx_dict[self.feature_type],
+ energies=energies,
+ e_range=self.e_range,
+ )
+ / self._get_n_moment(
+ n=2,
+ coxx=coxx_dict[self.feature_type],
+ energies=energies,
+ e_range=self.e_range,
+ )
+ ** 2
+ )
+ except KeyError:
+ raise ValueError(
+ "Please recheck fp_type requested argument.Possible options are bonding/antibonding/overall"
+ )
+
+ return coxx_center, coxx_width, coxx_skew, coxx_kurt
+
+[docs] def _get_coxx_center(
+ self,
+ coxx: npt.NDArray[np.floating],
+ energies: npt.NDArray[np.floating],
+ e_range: list[float],
+ ) -> float:
+ """
+ Get the bandwidth, defined as the first moment of the COXX.
+
+ Args:
+ coxx: COXX array
+ energies: energies corresponding COXX
+ e_range: range of energy to compute coxx center
+
+ Returns:
+ coxx center in eV
+ """
+ return self._get_n_moment(
+ n=1, coxx=coxx, energies=energies, e_range=e_range, center=False
+ )
+
+[docs] def _get_n_moment(
+ self,
+ n: float,
+ coxx: npt.NDArray[np.floating],
+ energies: npt.NDArray[np.floating],
+ e_range: list[float] | None,
+ center: bool = True,
+ ) -> float:
+ """
+ Get the nth moment of COXX.
+
+ Args:
+ n: The order for the moment
+ coxx: COXX array
+ energies: energies array
+ center: Take moments with respect to the COXX center
+ e_range: range of energy to compute nth moment
+
+ Returns:
+ COXX nth moment in eV
+ """
+ if e_range:
+ min_e = self.e_range[0]
+ max_e = self.e_range[1]
+ if min_e is None:
+ min_e = min(energies)
+ if max_e is None:
+ max_e = max(energies)
+ else:
+ min_e = min(energies)
+ max_e = max(energies)
+
+ tol = 1e-6
+
+ mask = (energies >= min_e - tol) & (energies <= max_e + tol)
+
+ coxx = coxx[mask]
+ energies = energies[mask]
+
+ if center:
+ coxx_center = self._get_coxx_center(
+ coxx=coxx, energies=energies, e_range=[min_e, max_e]
+ )
+ p = energies - coxx_center
+ else:
+ p = energies
+
+ return np.trapz(p**n * coxx, x=energies) / np.trapz(coxx, x=energies)
+
+[docs] def get_summarized_coxx_df(
+ self,
+ ids: str | None = None,
+ label_list: list[str] | None = None,
+ per_bond=True,
+ ) -> pd.DataFrame:
+ """
+ Get a pandas dataframe with COXX features.
+
+ Features consist of weighted ICOXX, effective interaction number and
+ moment features of COXX in the selected energy range.
+
+ Returns:
+ Returns a pandas dataframe with cohp/cobi/coop related features as per input file
+
+ """
+ if ids:
+ df = pd.DataFrame(index=[ids])
+ else:
+ ids = Path(self.path_to_coxxcar).parent.name
+ df = pd.DataFrame(index=[ids])
+
+ (
+ per_bnd_xx,
+ per_antibnd_xx,
+ w_icoxx,
+ ein_xx,
+ ) = self._calculate_wicoxx_ein()
+
+ if self.are_coops:
+ cc, cw, cs, ck = self._calc_moment_features(
+ label_list=label_list, per_bond=per_bond
+ )
+ df.loc[ids, "bnd_wICOOP"] = per_bnd_xx
+ df.loc[ids, "antibnd_wICOOP"] = per_antibnd_xx
+ df.loc[ids, "w_ICOOP"] = w_icoxx
+ df.loc[ids, "EIN_ICOOP"] = ein_xx
+ df.loc[ids, "center_COOP"] = cc
+ df.loc[ids, "width_COOP"] = cw
+ df.loc[ids, "skewness_COOP"] = cs
+ df.loc[ids, "kurtosis_COOP"] = ck
+ elif self.are_cobis:
+ cc, cw, cs, ck = self._calc_moment_features(
+ label_list=label_list, per_bond=per_bond
+ )
+ df.loc[ids, "bnd_wICOBI"] = per_bnd_xx
+ df.loc[ids, "antibnd_wICOBI"] = per_antibnd_xx
+ df.loc[ids, "w_ICOBI"] = w_icoxx
+ df.loc[ids, "EIN_ICOBI"] = ein_xx
+ df.loc[ids, "center_COBI"] = cc
+ df.loc[ids, "width_COBI"] = cw
+ df.loc[ids, "skewness_COBI"] = cs
+ df.loc[ids, "kurtosis_COBI"] = ck
+ else:
+ cc, cw, cs, ck = self._calc_moment_features(
+ label_list=label_list, per_bond=per_bond
+ )
+ df.loc[ids, "bnd_wICOHP"] = per_bnd_xx
+ df.loc[ids, "antibnd_wICOHP"] = per_antibnd_xx
+ df.loc[ids, "w_ICOHP"] = w_icoxx
+ df.loc[ids, "EIN_ICOHP"] = ein_xx
+ df.loc[ids, "center_COHP"] = cc
+ df.loc[ids, "width_COHP"] = cw
+ df.loc[ids, "skewness_COHP"] = cs
+ df.loc[ids, "kurtosis_COHP"] = ck
+
+ return df
+
+
+[docs]class FeaturizeCharges:
+ """
+ Class to compute ionicity from CHARGE.lobster data.
+
+ Args:
+ path_to_structure: path to POSCAR
+ path_to_charge : path to CHARGE.lobster (e.g., "CHARGE.lobster")
+ charge_type : set charge type used for computing ionicity. Possible options are "Mulliken" or "Loewdin"
+
+ Attributes:
+ get_df: pandas dataframe
+
+ """
+
+ def __init__(
+ self,
+ path_to_structure: str | Path,
+ path_to_charge: str | Path,
+ charge_type: str,
+ ):
+ """
+ Compute the Ionicity of the structure from CHARGE.lobster data.
+
+ Args:
+ path_to_structure: path to POSCAR
+ path_to_charge : path to CHARGE.lobster (e.g., "CHARGE.lobster")
+ charge_type : set charge type used for computing ionicity. Possible options are "Mulliken" or "Loewdin"
+
+ """
+ self.path_to_structure = path_to_structure
+ self.path_to_charge = path_to_charge
+ self.charge_type = charge_type
+
+[docs] def _calc_ionicity(self) -> float:
+ r"""
+ Calculate ionicity of the crystal structure based on quantum chemical charges.
+
+ \\[I_{\text{{Charges}}} = \frac{1}{{N_{\text{{Atoms}}}}}\\sum_{i=1}^{N_{\text{{Atoms}}} }
+ \\left(\frac{q_i}{v_{\text{{eff}},i}}\right)\\]
+
+ Returns:
+ Ionicity of the structure
+
+ """
+ chargeobj = Charge(filename=self.path_to_charge)
+ structure = Structure.from_file(self.path_to_structure)
+
+ if self.charge_type.lower() not in ["mulliken", "loewdin"]:
+ raise ValueError(
+ "Please check the requested charge_type. "
+ "Possible options are `Mulliken` or `Loewdin`"
+ )
+
+ ch_veff = []
+ tol = 1e-6
+ for i, j in enumerate(getattr(chargeobj, self.charge_type.capitalize())):
+ if (
+ j > tol
+ and not structure.species[i].is_transition_metal
+ and (
+ not structure.species[i].is_actinoid
+ and not structure.species[i].is_lanthanoid
+ )
+ ):
+ valence_elec = element(structure.species[i].value)
+ val = j / (valence_elec.nvalence() - 0)
+ ch_veff.append(val)
+
+ elif (
+ j < tol
+ and not structure.species[i].is_transition_metal
+ and (
+ not structure.species[i].is_actinoid
+ and not structure.species[i].is_lanthanoid
+ )
+ ):
+ valence_elec = element(structure.species[i].value)
+ val = j / (valence_elec.nvalence() - 8)
+ ch_veff.append(val)
+
+ elif j > tol and (
+ structure.species[i].is_transition_metal
+ or structure.species[i].is_actinoid
+ or structure.species[i].is_lanthanoid
+ ):
+ val = j / (structure.species[i].max_oxidation_state - 0)
+ ch_veff.append(val)
+
+ elif j < tol and (
+ structure.species[i].is_transition_metal
+ or structure.species[i].is_actinoid
+ or structure.species[i].is_lanthanoid
+ ):
+ val = j / (structure.species[i].min_oxidation_state - 8)
+ ch_veff.append(val)
+
+ return sum(ch_veff) / structure.num_sites
+
+[docs] def get_df(self, ids: str | None = None) -> pd.DataFrame:
+ """
+ Return a pandas dataframe with computed ionicity as columns.
+
+ Returns:
+ Returns a pandas dataframe with ionicity
+
+ """
+ if ids:
+ df = pd.DataFrame(index=[ids])
+ else:
+ ids = Path(self.path_to_charge).parent.name
+ df = pd.DataFrame(index=[ids])
+
+ if self.charge_type.lower() == "mulliken":
+ df.loc[ids, "Ionicity_Mull"] = self._calc_ionicity()
+ else:
+ df.loc[ids, "Ionicity_Loew"] = self._calc_ionicity()
+
+ return df
+
+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""Here classes and functions to plot Lobster outputs are provided."""
+
+from __future__ import annotations
+
+import typing
+from itertools import cycle
+from typing import TYPE_CHECKING, Any
+
+import matplotlib as mpl
+import numpy as np
+import plotly.graph_objs as go
+from matplotlib import pyplot as plt
+from pkg_resources import resource_filename
+from pymatgen.core import Structure
+from pymatgen.electronic_structure.cohp import Cohp, CompleteCohp, IcohpCollection
+from pymatgen.electronic_structure.core import Spin
+from pymatgen.electronic_structure.dos import LobsterCompleteDos
+from pymatgen.electronic_structure.plotter import CohpPlotter, DosPlotter
+
+if TYPE_CHECKING:
+ from lobsterpy.cohp.analyze import Analysis
+from lobsterpy.plotting import layout_dicts as ld
+
+base_style = resource_filename("lobsterpy.plotting", "lobsterpy_base.mplstyle")
+
+
+[docs]def get_style_list(
+ no_base_style: bool = False,
+ styles: list[str | dict[str, Any]] | None = None,
+ **kwargs,
+) -> list[str | dict[str, Any]]:
+ """Get *args for matplotlib.style from user input.
+
+ Args:
+ no_base_style: If true, do not include lobsterpy_base.mplstyle
+ styles: User-requested styles. These can be paths to mplstyle files,
+ the names of known (matplotlib-supplied) styles, or dicts of rcParam options.
+ **kwargs: matplotlib-style sheet keyword arguments
+
+ Remaining kwargs are collected as a dict and take the highest priority.
+ """
+ base = [] if no_base_style else [base_style]
+ if styles is None:
+ styles = []
+
+ return base + styles + [kwargs]
+
+
+[docs]class PlainCohpPlotter(CohpPlotter):
+ """
+ Modified Pymatgen CohpPlotter with styling removed.
+
+ This allows the styling to be manipulated more easily using matplotlib
+ style sheets.
+ """
+
+[docs] def get_plot(
+ self,
+ ax: mpl.axes.Axes | None = None,
+ xlim: tuple[float, float] | None = None,
+ ylim: tuple[float, float] | None = None,
+ plot_negative: bool | None = None,
+ integrated: bool = False,
+ invert_axes: bool = True,
+ sigma: float | None = None,
+ ):
+ """
+ Get a matplotlib plot showing the COHP.
+
+ Args:
+ ax: Existing Matplotlib Axes object to plot to.
+ xlim: Specifies the x-axis limits. Defaults to None for
+ automatic determination.
+ ylim: Specifies the y-axis limits. Defaults to None for
+ automatic determination.
+ plot_negative: It is common to plot -COHP(E) so that the
+ sign means the same for COOPs and COHPs. Defaults to None
+ for automatic determination: If are_coops is True, this
+ will be set to False, else it will be set to True.
+ integrated: Switch to plot ICOHPs. Defaults to False.
+ invert_axes: Put the energies onto the y-axis, which is
+ common in chemistry.
+ sigma: Standard deviation of Gaussian broadening applied to
+ population data. If this is unset (None) no broadening will be
+ added.
+
+ Returns:
+ A matplotlib object.
+
+ """
+ if self.are_coops and not self.are_cobis:
+ cohp_label = "COOP"
+ elif self.are_cobis and not self.are_coops:
+ cohp_label = "COBI"
+ elif self.are_cobis and self.are_coops:
+ raise ValueError(
+ "Plot data should not contain COBI and COOP data at same time"
+ )
+ else:
+ cohp_label = "COHP" + " (eV)"
+
+ if plot_negative is None:
+ plot_negative = (not self.are_coops) and (not self.are_cobis)
+
+ if integrated:
+ cohp_label = "I" + cohp_label
+
+ if plot_negative:
+ cohp_label = "$-$" + cohp_label
+
+ energy_label = "$E - E_f$ (eV)" if self.zero_at_efermi else "$E$ (eV)"
+
+ colors = mpl.rcParams["axes.prop_cycle"].by_key()["color"]
+ ncolors = len(colors)
+
+ if ax is None:
+ _, ax = plt.subplots()
+
+ allpts = []
+ keys = self._cohps.keys()
+ for i, key in enumerate(keys):
+ energies = self._cohps[key]["energies"]
+ populations = (
+ self._cohps[key]["COHP"]
+ if not integrated
+ else self._cohps[key]["ICOHP"]
+ )
+ for spin in [Spin.up, Spin.down]:
+ if spin in populations:
+ if invert_axes:
+ x = -populations[spin] if plot_negative else populations[spin]
+ y = energies
+ x = self._broaden(y, x, sigma=sigma)
+ else:
+ x = energies
+ y = -populations[spin] if plot_negative else populations[spin]
+ y = self._broaden(x, y, sigma=sigma)
+ allpts.extend(list(zip(x, y)))
+ if spin == Spin.up:
+ ax.plot(
+ x,
+ y,
+ color=colors[i % ncolors],
+ linestyle="-",
+ label=str(key),
+ )
+ else:
+ ax.plot(x, y, color=colors[i % ncolors], linestyle="--")
+
+ if xlim:
+ ax.set_xlim(xlim)
+ xlim = ax.get_xlim()
+ assert isinstance(xlim, tuple)
+
+ if ylim:
+ ax.set_ylim(ylim)
+ else:
+ relevanty = [p[1] for p in allpts if xlim[0] < p[0] < xlim[1]]
+ if relevanty:
+ plt.ylim((min(relevanty), max(relevanty)))
+
+ grid_like_line_kwargs = {
+ "color": mpl.rcParams["grid.color"],
+ "linewidth": mpl.rcParams["grid.linewidth"],
+ "linestyle": mpl.rcParams["grid.linestyle"],
+ "alpha": mpl.rcParams["grid.alpha"],
+ "zorder": 0,
+ }
+
+ if not invert_axes:
+ ax.axhline(**grid_like_line_kwargs)
+
+ if self.zero_at_efermi:
+ ax.axvline(**grid_like_line_kwargs)
+
+ else:
+ ax.axvline(self._cohps[key]["efermi"], **grid_like_line_kwargs)
+ else:
+ ax.axvline(**grid_like_line_kwargs)
+
+ if self.zero_at_efermi:
+ ax.axhline(**grid_like_line_kwargs)
+ else:
+ ax.axhline(self._cohps[key]["efermi"], **grid_like_line_kwargs)
+
+ if invert_axes:
+ plt.xlabel(cohp_label)
+ plt.ylabel(energy_label)
+ else:
+ plt.xlabel(energy_label)
+ plt.ylabel(cohp_label)
+
+ _ = ax.legend()
+
+ return plt
+
+[docs] @staticmethod
+ def _broaden(energies: np.ndarray, population: np.ndarray, sigma=None, cutoff=4.0):
+ """Broaden the spectrum with a given standard deviation.
+
+ The population is convolved with a normalised Gaussian kernel. This
+ requires the energy grid to be regularly-spaced.
+
+ Args:
+ energies: Regularly-spaced energy series
+ population: Population data for broadening
+ sigma: Standard deviation for Gaussian broadening. If sigma is None
+ then the input data is returned without any processing.
+ cutoff: Range cutoff for broadening kernel, as a multiple of sigma.
+
+ Return:
+ Broadened population
+
+ """
+ from scipy.signal import convolve
+ from scipy.stats import norm
+
+ if sigma is None:
+ return population
+
+ spacing = np.mean(np.diff(energies))
+ if not np.allclose(np.diff(energies), spacing, atol=1e-5):
+ raise ValueError(
+ "Energy grid is not regular, cannot broaden with "
+ "discrete convolution."
+ )
+
+ # Obtain symmetric mesh for broadening kernel, centered on zero
+ kernel_x = np.arange(0, cutoff * sigma + 0.5 * spacing, spacing)
+ kernel_x = np.concatenate([-kernel_x[-1:1:-1], kernel_x])
+
+ kernel = norm.pdf(kernel_x, scale=sigma)
+
+ return convolve(population, kernel, mode="same") / kernel.sum()
+
+
+[docs]class PlainDosPlotter(DosPlotter):
+ """
+ Modified Pymatgen DosPlotter with styling removed.
+
+ This allows the styling to be manipulated more easily using matplotlib
+ style sheets. It also adds additional functionalities to plotter
+ """
+
+ def __init__(
+ self, zero_at_efermi: bool = True, stack: bool = False, sigma=None, summed=False
+ ) -> None:
+ """
+ Generate COHP or COOP or COBI plots.
+
+ Args:
+ zero_at_efermi (bool): Whether to shift all Dos to have zero energy at the
+ fermi energy. Defaults to True.
+ stack (bool): Whether to plot the DOS as a stacked area graph
+ sigma (float): Specify a standard deviation for Gaussian smearing
+ the DOS for nicer looking plots. Defaults to None for no
+ smearing.
+ summed (bool): Whether to plot the summed DOS
+
+ """
+ self.zero_at_efermi = zero_at_efermi
+ self.stack = stack
+ self.sigma = sigma
+ self._norm_val = True
+ self._doses = {} # type: ignore
+ self.summed = summed
+
+[docs] def add_dos(self, label: str, dos: LobsterCompleteDos) -> None:
+ """Add a dos for plotting.
+
+ Args:
+ label: label for the DOS. Must be unique.
+ dos: LobsterCompleteDos object
+ """
+ if dos.norm_vol is None:
+ self._norm_val = False
+ energies = dos.energies
+ if self.summed:
+ if self.sigma:
+ smeared_densities = dos.get_smeared_densities(self.sigma)
+ if Spin.down in smeared_densities:
+ added_densities = (
+ smeared_densities[Spin.up] + smeared_densities[Spin.down]
+ )
+ densities = {Spin.up: added_densities}
+ else:
+ densities = smeared_densities
+ else:
+ densities = {Spin.up: dos.get_densities()}
+ else:
+ densities = (
+ dos.get_smeared_densities(self.sigma) if self.sigma else dos.densities
+ )
+
+ efermi = dos.efermi
+
+ self._doses[label] = {
+ "energies": energies,
+ "densities": densities,
+ "efermi": efermi,
+ }
+
+[docs] def add_site_orbital_dos(self, dos: LobsterCompleteDos, orbital, site_index):
+ """Add orbital dos at particular site.
+
+ Args:
+ dos: LobsterCompleteDos object
+ orbital: Orbitals name at the site. Must be unique.
+ site_index: site index in the structure
+ """
+ if dos.norm_vol is None:
+ self._norm_val = False
+ site = dos.structure.sites[site_index]
+ avail_orbs = list(dos.pdos[site])
+ if orbital not in avail_orbs and orbital != "all":
+ str_orbs = ", ".join(avail_orbs)
+ raise ValueError(
+ f"Requested orbital is not available for this site, "
+ f"available orbitals are {str_orbs}"
+ )
+
+ if orbital == "all":
+ for orb in avail_orbs:
+ dos_obj = dos.get_site_orbital_dos(site=site, orbital=orb)
+ label = site.species_string + str(site_index + 1) + f": {orb}"
+ energies = dos_obj.energies
+ if self.summed:
+ if self.sigma:
+ smeared_densities = dos_obj.get_smeared_densities(self.sigma)
+ if Spin.down in smeared_densities:
+ added_densities = (
+ smeared_densities[Spin.up]
+ + smeared_densities[Spin.down]
+ )
+ densities = {Spin.up: added_densities}
+ else:
+ densities = smeared_densities
+ else:
+ densities = {Spin.up: dos_obj.get_densities()}
+ else:
+ densities = (
+ dos_obj.get_smeared_densities(self.sigma)
+ if self.sigma
+ else dos_obj.densities
+ )
+
+ efermi = dos_obj.efermi
+
+ self._doses[label] = {
+ "energies": energies,
+ "densities": densities,
+ "efermi": efermi,
+ }
+ else:
+ dos_obj = dos.get_site_orbital_dos(site=site, orbital=orbital)
+ label = site.species_string + str(site_index + 1) + f": {orbital}"
+
+ energies = dos_obj.energies
+ if self.summed:
+ if self.sigma:
+ smeared_densities = dos_obj.get_smeared_densities(self.sigma)
+ if Spin.down in smeared_densities:
+ added_densities = (
+ smeared_densities[Spin.up] + smeared_densities[Spin.down]
+ )
+ densities = {Spin.up: added_densities}
+ else:
+ densities = smeared_densities
+ else:
+ densities = {Spin.up: dos_obj.get_densities()}
+ else:
+ densities = (
+ dos_obj.get_smeared_densities(self.sigma)
+ if self.sigma
+ else dos_obj.densities
+ )
+
+ efermi = dos_obj.efermi
+
+ self._doses[label] = {
+ "energies": energies,
+ "densities": densities,
+ "efermi": efermi,
+ }
+
+[docs] @typing.no_type_check
+ def get_plot(
+ self,
+ ax: mpl.axes.Axes | None = None,
+ xlim: tuple[float, float] | None = None,
+ ylim: tuple[float, float] | None = None,
+ invert_axes: bool = False,
+ beta_dashed: bool = False,
+ sigma: float | None = None,
+ ):
+ """
+ Get a matplotlib plot showing the COHP.
+
+ Args:
+ ax: Existing Matplotlib Axes object to plot to.
+ xlim: Specifies the x-axis limits. Defaults to None for
+ automatic determination.
+ ylim: Specifies the y-axis limits. Defaults to None for
+ automatic determination.
+ invert_axes: Put the energies onto the y-axis, which is
+ common in chemistry.
+ beta_dashed: Plots the beta spin channel with a dashed line. Defaults to False
+ sigma: Standard deviation of Gaussian broadening applied to
+ population data. If this is unset (None) no broadening will be
+ added.
+
+ Returns:
+ A matplotlib object.
+
+ """
+ ys = None
+ all_densities = []
+ all_energies = []
+
+ colors = mpl.rcParams["axes.prop_cycle"].by_key()["color"]
+ n_colors = len(colors)
+
+ if ax is None:
+ _, ax = plt.subplots()
+
+ # Note that this complicated processing of energies is to allow for
+ # stacked plots in matplotlib.
+ for dos in self._doses.values():
+ energies = dos["energies"]
+ densities = dos["densities"]
+ if not ys:
+ ys = {
+ Spin.up: np.zeros(energies.shape),
+ Spin.down: np.zeros(energies.shape),
+ }
+ new_dens = {}
+ for spin in [Spin.up, Spin.down]:
+ if spin in densities:
+ if self.stack:
+ ys[spin] += densities[spin]
+ new_dens[spin] = ys[spin].copy()
+ else:
+ new_dens[spin] = densities[spin]
+ all_energies.append(energies)
+ all_densities.append(new_dens)
+
+ keys = list(reversed(self._doses))
+ all_densities.reverse()
+ all_energies.reverse()
+ all_pts = []
+
+ for idx, key in enumerate(keys):
+ for spin in [Spin.up, Spin.down]:
+ if spin in all_densities[idx]:
+ energy = all_energies[idx]
+ densities = list(int(spin) * all_densities[idx][spin])
+ if invert_axes:
+ x = densities
+ y = energy
+ else:
+ x = energy
+ y = densities
+ all_pts.extend(list(zip(x, y)))
+ if self.stack:
+ ax.fill(x, y, color=colors[idx % n_colors], label=str(key))
+ elif spin == Spin.down and beta_dashed:
+ ax.plot(
+ x,
+ y,
+ color=colors[idx % n_colors],
+ label=str(key),
+ linestyle="--",
+ linewidth=3,
+ )
+ else:
+ ax.plot(
+ x,
+ y,
+ color=colors[idx % n_colors],
+ label=str(key),
+ linewidth=3,
+ )
+
+ if xlim:
+ ax.set_xlim(xlim)
+ if ylim:
+ ax.set_ylim(ylim)
+ elif not invert_axes:
+ xlim = ax.get_xlim()
+ relevant_y = [p[1] for p in all_pts if xlim[0] < p[0] < xlim[1]]
+ ax.set_ylim((min(relevant_y), max(relevant_y)))
+ if not xlim and invert_axes:
+ ylim = ax.get_ylim()
+ relevant_y = [p[0] for p in all_pts if ylim[0] < p[1] < ylim[1]]
+ ax.set_xlim((min(relevant_y), max(relevant_y)))
+
+ if self.zero_at_efermi:
+ xlim = ax.get_xlim()
+ ylim = ax.get_ylim()
+ ax.plot(xlim, [0, 0], "k--", linewidth=2) if invert_axes else ax.plot(
+ [0, 0], ylim, "k--", linewidth=2
+ )
+
+ if invert_axes:
+ ax.set_ylabel("Energies (eV)")
+ ax.set_xlabel(
+ f"Density of states (states/eV{'/ų' if self._norm_val else ''})"
+ )
+ ax.axvline(x=0, color="k", linestyle="--", linewidth=2)
+ else:
+ ax.set_xlabel("Energies (eV)")
+ if self._norm_val:
+ ax.set_ylabel("Density of states (states/eV/ų)")
+ else:
+ ax.set_ylabel("Density of states (states/eV)")
+ ax.axhline(y=0, color="k", linestyle="--", linewidth=2)
+
+ # Remove duplicate labels with a dictionary
+ handles, labels = ax.get_legend_handles_labels()
+ label_dict = dict(zip(labels, handles))
+ ax.legend(label_dict.values(), label_dict.keys())
+ legend_text = (
+ ax.get_legend().get_texts()
+ ) # all the text.Text instance in the legend
+ plt.setp(legend_text, fontsize=30)
+ plt.tight_layout()
+ _ = ax.legend()
+
+ return plt
+
+
+[docs]class InteractiveCohpPlotter(CohpPlotter):
+ """Interactive COHP, COBI or COOP plotter to view all relevant bonds in one figure."""
+
+ COLOR_PALETTE = [
+ "#e41a1c",
+ "#377eb8",
+ "#4daf4a",
+ "#984ea3",
+ "#ff7f00",
+ "#ffff33",
+ "#a65628",
+ "#f781bf",
+ "#999999",
+ ]
+
+[docs] def add_all_relevant_cohps(
+ self,
+ analyse: Analysis,
+ suffix: str = "",
+ label_resolved: bool = True,
+ orbital_resolved: bool = False,
+ ) -> None:
+ """
+ Add all relevant COHPs from lobsterpy analyse object.
+
+ Args:
+ analyse: Analyse object from lobsterpy.
+ suffix: Optional addition to LOBSTER label to avoid key conflicts when plotting multiple
+ calcs or just for additional legend information.
+ label_resolved: bool indicating to obtain label resolved interactive plots for relevant bonds.
+ If false, will return summed cohp curves of unique relevant bonds.
+ orbital_resolved: bool indicating to include orbital resolved interactive cohps for relevant bonds.
+ """
+ complete_cohp = analyse.chemenv.completecohp
+
+ plot_data = analyse.get_site_bond_resolved_labels()
+
+ if "All" not in self._cohps:
+ self._cohps["All"] = {}
+
+ # iterate and extract the data to be plotted from cohp objects
+ for bond_key, labels in plot_data.items():
+ count = len(labels)
+ label_with_count = self._insert_number_of_bonds_in_label(
+ label=bond_key, character=":", number_of_bonds=count
+ )
+ # get summed cohps from the relevant bond label at the site
+ cohp = complete_cohp.get_summed_cohp_by_label_list(label_list=labels)
+ energies = (
+ cohp.energies - cohp.efermi if self.zero_at_efermi else cohp.energies
+ )
+ drop_down_key = plot_legend = (
+ label_with_count + suffix
+ ) # label for the dropdown menu
+ self._cohps["All"].update(
+ {
+ plot_legend: {
+ "energies": energies,
+ "COHP": cohp.get_cohp(),
+ "ICOHP": cohp.get_icohp(),
+ "efermi": cohp.efermi,
+ }
+ }
+ )
+ if len(plot_data) > 1:
+ if drop_down_key not in self._cohps:
+ self._cohps[drop_down_key] = {}
+ self._cohps[drop_down_key].update(
+ {
+ plot_legend: {
+ "energies": energies,
+ "COHP": cohp.get_cohp(),
+ "ICOHP": cohp.get_icohp(),
+ "efermi": cohp.efermi,
+ }
+ }
+ )
+
+ # Add cohp data for each relevant bond label iteratively
+ if label_resolved and not orbital_resolved:
+ if label_with_count + suffix not in self._cohps:
+ self._cohps[label_with_count + suffix] = {}
+ # Get cohp data for each relevant bond label at the site
+ for label in labels:
+ cohp = complete_cohp.get_cohp_by_label(label)
+ energies = (
+ cohp.energies - cohp.efermi
+ if self.zero_at_efermi
+ else cohp.energies
+ )
+ drop_down_key = (
+ label_with_count + suffix
+ ) # label for the dropdown menu
+ plot_legend = self._get_plot_label_for_label_resolved(
+ structure=analyse.structure,
+ label_list=[label],
+ complete_cohp=complete_cohp,
+ orb_list=[],
+ label_resolved=True,
+ orbital_resolved=False,
+ )
+ if len(plot_data) > 1:
+ self._cohps[drop_down_key].update(
+ {
+ plot_legend: {
+ "energies": energies,
+ "COHP": cohp.get_cohp(),
+ "ICOHP": cohp.get_icohp(),
+ "efermi": cohp.efermi,
+ }
+ }
+ )
+
+ plot_legend_here = plot_legend + suffix
+ self._cohps["All"].update(
+ {
+ plot_legend_here: {
+ "energies": energies,
+ "COHP": cohp.get_cohp(),
+ "ICOHP": cohp.get_icohp(),
+ "efermi": cohp.efermi,
+ }
+ }
+ )
+ # Adds cohp data for each relevant orbitals and bond label iteratively
+ if orbital_resolved and label_resolved:
+ if label_with_count + suffix not in self._cohps:
+ self._cohps[label_with_count + suffix] = {}
+ # get relevant orbitals associated with each bond label
+ plot_data_orb = analyse.get_site_orbital_resolved_labels()
+ drop_down_key = label_with_count + suffix # label for the dropdown menu
+ key_val = plot_data_orb[bond_key]
+ # get cohp data for each orbital and associated bond labels iteratively
+ for orb, val in key_val.items():
+ for lab in val:
+ cohp_orb = (
+ complete_cohp.get_summed_cohp_by_label_and_orbital_list(
+ label_list=[lab], orbital_list=[orb]
+ )
+ )
+
+ energies = (
+ cohp_orb.energies - cohp_orb.efermi
+ if self.zero_at_efermi
+ else cohp_orb.energies
+ )
+ # plot legends will contain species and orbital along with relevant bond label
+ plot_legend = self._get_plot_label_for_label_resolved(
+ structure=analyse.structure,
+ label_list=[lab],
+ complete_cohp=complete_cohp,
+ orb_list=[orb],
+ orbital_resolved=True,
+ label_resolved=True,
+ )
+ if len(plot_data) > 1:
+ self._cohps[drop_down_key].update(
+ {
+ plot_legend: {
+ "energies": energies,
+ "COHP": cohp_orb.get_cohp(),
+ "ICOHP": cohp_orb.get_icohp(),
+ "efermi": cohp_orb.efermi,
+ }
+ }
+ )
+
+ plot_legend_here = plot_legend + suffix
+
+ self._cohps["All"].update(
+ {
+ plot_legend_here: {
+ "energies": energies,
+ "COHP": cohp_orb.get_cohp(),
+ "ICOHP": cohp_orb.get_icohp(),
+ "efermi": cohp_orb.efermi,
+ }
+ }
+ )
+ # Adds summed cohp data for each relevant orbitals
+ if orbital_resolved and not label_resolved:
+ if label_with_count + suffix not in self._cohps:
+ self._cohps[label_with_count + suffix] = {}
+ # get relevant orbitals associated with each bond label
+ plot_data_orb = analyse.get_site_orbital_resolved_labels()
+ drop_down_key = label_with_count + suffix # label for the dropdown menu
+ key_val = plot_data_orb[bond_key]
+ # get summed cohp data for each relevant orbital
+ for orb, val in key_val.items():
+ cohp_orb = complete_cohp.get_summed_cohp_by_label_and_orbital_list(
+ label_list=val, orbital_list=[orb] * len(val)
+ )
+
+ energies = (
+ cohp_orb.energies - cohp_orb.efermi
+ if self.zero_at_efermi
+ else cohp_orb.energies
+ )
+ # plot legends will contain species and orbital along with number of bonds at the site
+ plot_legend = self._get_plot_label_for_label_resolved(
+ structure=analyse.structure,
+ label_list=val,
+ complete_cohp=complete_cohp,
+ orb_list=[orb],
+ orbital_resolved=True,
+ label_resolved=False,
+ )
+
+ if len(plot_data) > 1:
+ self._cohps[drop_down_key].update(
+ {
+ plot_legend: {
+ "energies": energies,
+ "COHP": cohp_orb.get_cohp(),
+ "ICOHP": cohp_orb.get_icohp(),
+ "efermi": cohp_orb.efermi,
+ }
+ }
+ )
+
+ plot_legend_here = plot_legend + suffix
+
+ self._cohps["All"].update(
+ {
+ plot_legend_here: {
+ "energies": energies,
+ "COHP": cohp_orb.get_cohp(),
+ "ICOHP": cohp_orb.get_icohp(),
+ "efermi": cohp_orb.efermi,
+ }
+ }
+ )
+
+[docs] def add_cohps_by_lobster_label(
+ self, analyse: Analysis, label_list: list, suffix: str = ""
+ ):
+ """
+ Add COHPs explicitly specified in label list.
+
+ Args:
+ analyse: Analyse object from lobsterpy.
+ label_list: List of COHP labels as from LOBSTER.
+ suffix: Optional addition to LOBSTER label to avoid key
+ conflicts when plotting multiple calcs or just for additional legend information.
+ """
+ complete_cohp = analyse.chemenv.completecohp
+
+ if "All" not in self._cohps:
+ self._cohps["All"] = {}
+
+ for label in label_list:
+ atom1 = complete_cohp.bonds[label]["sites"][0].species_string
+ atom2 = complete_cohp.bonds[label]["sites"][1].species_string
+ sorted_label = sorted([atom1, atom2])
+ cohp = complete_cohp.get_cohp_by_label(label)
+ energies = (
+ cohp.energies - cohp.efermi if self.zero_at_efermi else cohp.energies
+ )
+ key = sorted_label[0] + "-" + sorted_label[1] + ": " + label + suffix
+ self._cohps["All"].update(
+ {
+ key: {
+ "energies": energies,
+ "COHP": cohp.get_cohp(),
+ "ICOHP": cohp.get_icohp(),
+ "efermi": cohp.efermi,
+ }
+ }
+ )
+
+[docs] def add_cohps_from_plot_data(self, plot_data_dict: dict, suffix: str = ""):
+ """
+ Add all relevant COHPs for specified bond type from lobster lightweight json.gz file.
+
+ Args:
+ plot_data_dict: Lobsterpy plot data dict
+ suffix: Optional addition to LOBSTER label to avoid key
+ conflicts when plotting multiple calcs or just for additional legend information.
+ """
+ # convert to cohp objects
+ plot_data_dict = plot_data_dict.copy()
+ for key, cohps in plot_data_dict.items():
+ if isinstance(cohps, Cohp):
+ plot_data_dict.update({key: cohps})
+ else:
+ try:
+ cohps = Cohp.from_dict(cohps)
+ plot_data_dict.update({key: cohps})
+ except TypeError:
+ raise ValueError(
+ "The data provided could not be converted to cohp object.Please recheck the input data"
+ )
+
+ if "All" not in self._cohps:
+ self._cohps["All"] = {}
+
+ for bond_key, cohps in plot_data_dict.items():
+ energies = (
+ cohps.energies - cohps.efermi if self.zero_at_efermi else cohps.energies
+ )
+ key = bond_key + suffix
+ self._cohps["All"].update(
+ {
+ key: {
+ "energies": energies,
+ "COHP": cohps.get_cohp(),
+ "ICOHP": cohps.get_icohp(),
+ "efermi": cohps.efermi,
+ }
+ }
+ )
+
+[docs] def get_plot(
+ self,
+ xlim=None,
+ rangeslider=False,
+ ylim=None,
+ plot_negative=None,
+ integrated=False,
+ invert_axes=True,
+ sigma=None,
+ colors=None,
+ ):
+ """
+ Get an interactive plotly figure showing the COHPs.
+
+ Args:
+ xlim: Specifies the x-axis limits. Defaults to None for
+ automatic determination.
+ rangeslider: Adds a plotly.graph_objs.layout.xaxis.Rangeslider
+ object to figure to allow easy manipulation of x-axis..
+ ylim: Specifies the y-axis limits. Defaults to None for
+ automatic determination.
+ plot_negative: It is common to plot -COHP(E) so that the
+ sign means the same for COOPs and COHPs. Defaults to None
+ for automatic determination: If are_coops is True, this
+ will be set to False, else it will be set to True.
+ integrated: Switch to plot ICOHPs. Defaults to False.
+ invert_axes: Put the energies onto the y-axis, which is
+ common in chemistry.
+ sigma: Standard deviation of Gaussian broadening applied to
+ population data. If this is unset (None) no broadening will be added.
+ colors: list of hex color codes to be used in plot
+
+ Returns:
+ A plotly.graph_objects.Figure object.
+
+ """
+ if self.are_coops and not self.are_cobis:
+ cohp_label = "COOP"
+ elif self.are_cobis and not self.are_coops:
+ cohp_label = "COBI"
+ elif self.are_cobis and self.are_coops:
+ raise ValueError(
+ "Plot data should not contain COBI and COOP data at same time"
+ )
+ else:
+ cohp_label = "COHP" + " (eV)"
+
+ if plot_negative is None:
+ plot_negative = (not self.are_coops) and (not self.are_cobis)
+
+ if integrated:
+ cohp_label = "I" + cohp_label
+
+ if plot_negative:
+ cohp_label = "\u2212" + cohp_label
+
+ energy_label = (
+ "$E - E_f \\text{ (eV)}$" if self.zero_at_efermi else "$E \\text{ (eV)}$"
+ )
+
+ # Setting up repeating color scheme (same as for matplotlib plots in .mplstyle)
+ palette = InteractiveCohpPlotter.COLOR_PALETTE if colors is None else colors
+ pal_iter = cycle(palette)
+
+ traces = {}
+ for k, v in self._cohps.items():
+ traces.update({k: []})
+ for label in v:
+ population_key = v[label]["ICOHP"] if integrated else v[label]["COHP"]
+ band_color = next(pal_iter)
+ for spin in [Spin.up, Spin.down]:
+ if spin in population_key:
+ population = (
+ [-i for i in population_key[spin]]
+ if plot_negative
+ else population_key[spin]
+ )
+ if invert_axes:
+ x = population
+ y = v[label]["energies"]
+ x = PlainCohpPlotter._broaden(y, x, sigma=sigma)
+ else:
+ x = v[label]["energies"]
+ y = population
+ y = PlainCohpPlotter._broaden(x, y, sigma=sigma)
+ if spin == Spin.up:
+ trace = go.Scatter(x=x, y=y, name=label)
+ trace.update(ld.spin_up_trace_style_dict)
+ else:
+ trace = go.Scatter(x=x, y=y, name="")
+ trace.update(ld.spin_down_trace_style_dict)
+ trace.update(line={"color": band_color})
+ traces[k].append(trace)
+
+ energy_axis = (
+ go.layout.YAxis(title=energy_label)
+ if invert_axes
+ else go.layout.XAxis(
+ title=energy_label, rangeslider={"visible": rangeslider}
+ )
+ )
+ energy_axis.update(ld.energy_axis_style_dict)
+ cohp_axis = (
+ go.layout.XAxis(title=cohp_label, rangeslider={"visible": rangeslider})
+ if invert_axes
+ else go.layout.YAxis(title=cohp_label)
+ )
+ cohp_axis.update(ld.cohp_axis_style_dict)
+
+ layout = (
+ go.Layout(xaxis=cohp_axis, yaxis=energy_axis)
+ if invert_axes
+ else go.Layout(xaxis=energy_axis, yaxis=cohp_axis)
+ )
+
+ # Create figure object
+ fig = go.Figure(layout=layout)
+ fig.update_layout(ld.layout_dict)
+ fig.update_layout(legend=ld.legend_style_dict)
+
+ # Add all traces to figure
+ for val_trace in traces.values():
+ for trace in val_trace:
+ fig.add_trace(trace)
+
+ # Update layout with dropdown menu if label resolved plot
+ if len(traces) > 2:
+ # Update visibility of traces
+ for i, _ in enumerate(fig.data):
+ if i <= len(traces["All"]) - 1:
+ pass
+ else:
+ fig.data[i].visible = False
+ # Add dropdown buttons
+ fig.update_layout(
+ updatemenus=[
+ {
+ "buttons": [
+ {
+ "args": [
+ {
+ "visible": [
+ selected_group == group
+ for group, val_trace in traces.items()
+ for _trace in val_trace
+ ]
+ }
+ ],
+ "label": selected_group,
+ "method": "update",
+ }
+ for selected_group in traces
+ ],
+ "direction": "down",
+ "showactive": True,
+ "active": 0,
+ "x": 0.5,
+ "y": 1.15,
+ "bgcolor": "rgba(255,255,255,0.8)",
+ "bordercolor": "rgba(0,0,0,0.2)",
+ "xanchor": "center",
+ "yanchor": "top",
+ "font": {"family": "Arial", "color": "#444444", "size": 18},
+ }
+ ]
+ )
+
+ if xlim:
+ fig.update_xaxes(range=xlim)
+ if ylim:
+ fig.update_yaxes(range=ylim)
+
+ fig.update_yaxes(automargin=True)
+
+ return fig
+
+[docs] @staticmethod
+ def _insert_number_of_bonds_in_label(
+ label: str, character: str, number_of_bonds: int
+ ) -> str:
+ """
+ Add number of bonds to bond label.
+
+ For example : for input label 'Ba1: Ba-Ti', character ':', number_of_bonds: 3,
+ Will return 'Ba1: 3 x Ba-Ti'
+
+ Args:
+ label: bond label to which number of bonds needs to be inserted
+ character: string character where number of bonds needs to be inserted
+ number_of_bonds: number of bonds corresponding to the label
+
+ Returns:
+ bond label with number of bonds inserted
+ """
+ return label.replace(character, f"{character} {number_of_bonds} x", 1)
+
+[docs] @staticmethod
+ def _get_plot_label_for_label_resolved(
+ structure: Structure,
+ label_list: list,
+ complete_cohp: CompleteCohp,
+ orb_list: list,
+ label_resolved: bool = False,
+ orbital_resolved: bool = False,
+ ) -> str:
+ """
+ Generate plot labels for both orbital and label-resolved plots.
+
+ For example for NaCl structure, label:21, orb: 3s-3s
+ Will return '21: Cl2-Na1 (2.85 Å)' if label_resolved is True and orbital_resolved is False
+ Will return '21: Cl2(3s)-Na1(3s) (2.85 Å)' If label and orbital resolved True
+ Will return 'Cl(3s)-Na(3s) (2.85 Å)' if orbital_resolved is True and label_resolved is False
+
+ Args:
+ structure: pymatgen structure object
+ label_list: bond label to which number of bonds needs to be inserted
+ complete_cohp: complete cohp object
+ orb_list: relevant orbital
+ label_resolved: specifies type of label to be returned is for label_resolved case
+ orbital_resolved: specifies type of label to be returned is for orbital_resolved case
+
+ Returns:
+ plot label string
+ """
+ if label_resolved and not orbital_resolved:
+ atom_pairs = []
+ for site in complete_cohp.bonds[label_list[0]]["sites"]:
+ atom = f"{site.species_string}{structure.sites.index(site) + 1!s}"
+ atom_pairs.append(atom)
+ atom_pair_str = "-".join(atom_pairs)
+ bond_length = round(complete_cohp.bonds[label_list[0]]["length"], 2)
+ plot_label = f"{label_list[0]}: {atom_pair_str} ({bond_length} \u00c5)"
+
+ elif not label_resolved and orbital_resolved:
+ orb_atom_pairs = []
+ orb_pair = orb_list[0].split("-")
+ for site, site_orb in zip(
+ complete_cohp.bonds[label_list[0]]["sites"], orb_pair
+ ):
+ atom_orb = f"{site.species_string} ({site_orb})"
+ orb_atom_pairs.append(atom_orb)
+ orb_atom_pairs_str = "-".join(orb_atom_pairs)
+ bond_length = round(complete_cohp.bonds[label_list[0]]["length"], 2)
+ plot_label = (
+ f"{len(label_list)}x {orb_atom_pairs_str} ({bond_length} \u00c5)"
+ )
+ else:
+ orb_atom_pairs = []
+ orb_pair = orb_list[0].split("-")
+ for site, site_orb in zip(
+ complete_cohp.bonds[label_list[0]]["sites"], orb_pair
+ ):
+ atom_orb = f"{site.species_string}{structure.sites.index(site) + 1!s} ({site_orb})"
+ orb_atom_pairs.append(atom_orb)
+ orb_atom_pairs_str = "-".join(orb_atom_pairs)
+ bond_length = round(complete_cohp.bonds[label_list[0]]["length"], 2)
+ plot_label = (
+ f"{label_list[0]}: {orb_atom_pairs_str} ({bond_length} \u00c5)"
+ )
+
+ return plot_label
+
+
+[docs]class IcohpDistancePlotter:
+ """Plotter to generate ICOHP or ICOBI or ICOOP vs bond lengths plots."""
+
+ def __init__(self, are_coops: bool = False, are_cobis: bool = False):
+ """
+ Plot ICOHPs or ICOBI or ICOOP vs bond lengths.
+
+ Args:
+ are_coops: Switch to indicate that these are ICOOPs, not ICOHPs.
+ Defaults to False for ICOHPs.
+ are_cobis: Switch to indicate that these are ICOBIs, not ICOHPs/COOPs.
+ Defaults to False for ICOHPs.
+ """
+ self.are_coops = are_coops
+ self.are_cobis = are_cobis
+ self._icohps = {} # type: ignore
+
+[docs] def add_icohps(self, label, icohpcollection: IcohpCollection):
+ """
+ Add ICOHPs or ICOBIs or ICOOPS for plotting.
+
+ Args:
+ label: Label for the ICOHPs. Must be unique.
+ icohpcollection: IcohpCollection object.
+
+ """
+ icohps = []
+ bond_len = []
+ atom_pairs = []
+ orb_data = {} # type: ignore
+ for indx, bond_label in enumerate(icohpcollection._list_labels):
+ orb_data.update({bond_label: {}})
+ for k, v in icohpcollection._list_orb_icohp[indx].items():
+ orb_data[bond_label].update({k: sum(v["icohp"].values())})
+ icohps.append(sum(icohpcollection._list_icohp[indx].values()))
+ bond_len.append(icohpcollection._list_length[indx])
+ atom1 = icohpcollection._list_atom1[indx]
+ atom2 = icohpcollection._list_atom2[indx]
+ atom_pairs.append(atom1 + "-" + atom2)
+
+ self._icohps[label] = {
+ "atom_pairs": atom_pairs,
+ "bond_labels": icohpcollection._list_labels,
+ "icohps": icohps,
+ "bond_lengths": bond_len,
+ "orb_data": orb_data,
+ }
+
+[docs] def get_plot(
+ self,
+ ax: mpl.axes.Axes | None = None,
+ marker_size: float = 50,
+ marker_style: str = "o",
+ xlim: tuple[float, float] | None = None,
+ ylim: tuple[float, float] | None = None,
+ plot_negative: bool = True,
+ ):
+ """
+ Get a matplotlib plot showing the COHP or COBI or COOP with respect to bond lengths.
+
+ Args:
+ ax: Existing Matplotlib Axes object to plot to.
+ marker_size: sets the size of markers in scatter plots
+ marker_style: sets type of marker used in plot
+ xlim: Specifies the x-axis limits. Defaults to None for
+ automatic determination.
+ ylim: Specifies the y-axis limits. Defaults to None for
+ automatic determination.
+ plot_negative: Will plot -1*ICOHPs. Works only for ICOHPs
+
+ Returns:
+ A matplotlib object.
+ """
+ if self.are_coops and not self.are_cobis:
+ cohp_label = "ICOOP"
+ elif self.are_cobis and not self.are_coops:
+ cohp_label = "ICOBI"
+ elif self.are_cobis and self.are_coops:
+ raise ValueError(
+ "Plot data should not contain ICOBI and ICOOP data at same time"
+ )
+ else:
+ cohp_label = "ICOHP (eV)"
+
+ if ax is None:
+ _, ax = plt.subplots()
+
+ if xlim:
+ ax.set_xlim(xlim)
+ if ylim:
+ ax.set_ylim(ylim)
+
+ ax.set_ylabel(cohp_label)
+ ax.set_xlabel("Bond lengths (\u00c5)")
+
+ for label, data in self._icohps.items():
+ x = data["bond_lengths"]
+ if plot_negative and cohp_label == "ICOHP (eV)":
+ ax.set_ylabel("$-$" + cohp_label)
+ y = [-1 * icohp for icohp in data["icohps"]]
+ else:
+ y = data["icohps"]
+
+ ax.scatter(x, y, s=marker_size, marker=marker_style, label=label)
+
+ return plt
+
+# Copyright (c) lobsterpy development team
+# Distributed under the terms of a BSD 3-Clause "New" or "Revised" License
+
+"""This module provides a class for generating graph objects using lobsterpy data."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+if TYPE_CHECKING:
+ from pathlib import Path
+
+from pymatgen.core.structure import Structure
+from pymatgen.io.lobster.lobsterenv import LobsterNeighbors
+from pymatgen.io.lobster.outputs import Charge
+
+from lobsterpy.cohp.analyze import Analysis
+
+
+[docs]class LobsterGraph:
+ """
+ Class to generate structure graph objects with bonding data from Lobster.
+
+ Attributes:
+ sg: return structure_graph object
+ """
+
+ def __init__(
+ self,
+ path_to_poscar: str | Path,
+ path_to_charge: str | Path,
+ path_to_cohpcar: str | Path,
+ path_to_icohplist: str | Path,
+ path_to_madelung: str | Path,
+ add_additional_data_sg=True,
+ path_to_icooplist: str | Path | None = None,
+ path_to_icobilist: str | Path | None = None,
+ which_bonds: str = "all",
+ cutoff_icohp: float = 0.10,
+ start: float | None = None,
+ ):
+ """
+ Return a structure graph objects with bonding information from Lobster data.
+
+ Mode of automatic bonding analysis can be “cation-anion” or “all” bonds. The strongest bond is
+ determined based on the ICOHPs. The coordination environments are determined based on
+ cutoff_icohp *ICOHPs values. If the path of ICOBILIST (ICOOPLIST) is provided, the ICOBI (ICOOP)
+ values corresponding to relevant bond labels obtained from the ICOHPLIST are also added as edge properties
+ to the structure graph objects. The Mulliken and Loewdin charges are added as node properties to
+ the structure graph objects.
+
+
+ Args:
+ path_to_poscar: path to POSCAR (e.g., "POSCAR")
+ path_to_charge: path to CHARGE.lobster (e.g., "CHARGE.lobster")
+ path_to_cohpcar: path to COHPCAR.lobster (e.g., "COHPCAR.lobster")
+ path_to_icohplist: path to ICOHPLIST.lobster (e.g., "ICOHPLIST.lobster")
+ path_to_icooplist: path to ICOOPLIST.lobster (e.g., "ICOOPLIST.lobster")
+ path_to_icobilist: path to ICOBILIST.lobster (e.g., "ICOBILIST.lobster")
+ path_to_madelung: path to MadelungEnergies.lobster (e.g., "MadelungEnergies.lobster")
+ cutoff_icohp : only bonds that are stronger than cutoff_icohp*strongest ICOHP will be considered
+ add_additional_data_sg: (bool) if True will add the information from ICOOPLIST.lobster
+ and ICOBILIST.lobster based on ICOHPLIST.lobster relevant bond
+ which_bonds: selects which kind of bonds are analyzed. "all" is the default
+ start: start energy for bonding antibonding percent integration
+
+ """
+ if add_additional_data_sg:
+ self.add_additional_data_sg = add_additional_data_sg
+ if path_to_icooplist is not None and path_to_icobilist is not None:
+ self.path_to_icooplist = path_to_icooplist
+ self.path_to_icobilist = path_to_icobilist
+ else:
+ raise ValueError(
+ "add_additional_data_sg is set to True."
+ "Please provide path_to_icooplist and path_to_icobilist"
+ )
+ else:
+ self.add_additional_data_sg = add_additional_data_sg
+
+ self.path_to_poscar = path_to_poscar
+ self.path_to_charge = path_to_charge
+ self.path_to_cohpcar = path_to_cohpcar
+ self.path_to_icohplist = path_to_icohplist
+ self.path_to_madelung = path_to_madelung
+ self.which_bonds = which_bonds
+ self.cutoff_icohp = cutoff_icohp
+
+ if self.which_bonds == "all":
+ self.additional_condition = 0
+ elif self.which_bonds == "cation-anion":
+ self.additional_condition = 1
+ else:
+ raise ValueError(
+ "Only accepted values are 'all' and 'cation-anion'."
+ "Please check the input parameters of which_bonds arg"
+ )
+ self.start = start
+
+ self.sg = self.get_decorated_sg()
+
+[docs] def get_decorated_sg(self):
+ """
+ Create a graph object decorated with bonding data from LobsterPy.
+
+ Returns:
+ structure graph object
+ """
+ if self.add_additional_data_sg:
+ chemenvlobster = LobsterNeighbors(
+ are_coops=False,
+ filename_ICOHP=self.path_to_icohplist,
+ perc_strength_ICOHP=self.cutoff_icohp,
+ structure=Structure.from_file(self.path_to_poscar),
+ additional_condition=self.additional_condition,
+ filename_CHARGE=self.path_to_charge,
+ add_additional_data_sg=self.add_additional_data_sg,
+ filename_blist_sg1=self.path_to_icobilist,
+ id_blist_sg1="ICOBI",
+ filename_blist_sg2=self.path_to_icooplist,
+ id_blist_sg2="ICOOP",
+ valences_from_charges=True,
+ adapt_extremum_to_add_cond=True,
+ )
+
+ else:
+ chemenvlobster = LobsterNeighbors(
+ are_coops=False,
+ filename_ICOHP=self.path_to_icohplist,
+ perc_strength_ICOHP=self.cutoff_icohp,
+ structure=Structure.from_file(self.path_to_poscar),
+ additional_condition=self.additional_condition,
+ filename_CHARGE=self.path_to_charge,
+ add_additional_data_sg=self.add_additional_data_sg,
+ valences_from_charges=True,
+ adapt_extremum_to_add_cond=True,
+ )
+
+ # Adds Mulliken and Löwdin charges as site properties to structure object (node properties)
+ decorated_structure = Charge(self.path_to_charge).get_structure_with_charges(
+ self.path_to_poscar
+ )
+
+ # Create the structure graph object decorated with site and edge properties based on ICOHP/ICOBI/ICOOP data
+ lobster_env = chemenvlobster.get_bonded_structure(
+ structure=decorated_structure, decorate=True, edge_properties=True
+ )
+
+ # Initialize automating bonding analysis from Lobsterpy based on ICOHP
+ analyze = Analysis(
+ path_to_charge=self.path_to_charge,
+ path_to_cohpcar=self.path_to_cohpcar,
+ path_to_poscar=self.path_to_poscar,
+ path_to_icohplist=self.path_to_icohplist,
+ path_to_madelung=self.path_to_madelung,
+ which_bonds=self.which_bonds,
+ cutoff_icohp=self.cutoff_icohp,
+ )
+
+ # Store the summarized dictionary object containing bonding information
+ cba = analyze.condensed_bonding_analysis
+
+ # Iterate over sites in the dictionary
+ for cba_data in cba["sites"].values():
+ for _, node_data in lobster_env.graph.nodes.data():
+ # Check if ions are same and add its corresponding environment in node properties
+ if cba_data["ion"] == node_data["specie"]:
+ node_data["properties"].update({"env": cba_data["env"]})
+
+ for (
+ edge_prop
+ ) in lobster_env.graph.edges.data(): # Iterate over structure graph edges
+ _ab, ab_p, _b, b_p = analyze._integrate_antbdstates_below_efermi(
+ cohp=analyze.chemenv.completecohp.get_cohp_by_label(
+ edge_prop[2]["bond_label"]
+ ),
+ start=self.start,
+ ) # Compute bonding- antibonding percentages for each bond in structure graph object
+ edge_prop[2][
+ "ICOHP_bonding_perc"
+ ] = b_p # Store bonding percentage in edge of graph object
+ edge_prop[2][
+ "ICOHP_antibonding_perc"
+ ] = ab_p # Store anti-bonding percentage in edge graph object
+
+ return lobster_env
+