Merge pull request #137 from naik-aakash/icohp_vs_bond_lengths

ICOHP vs bond length plotter

.github/workflows/python-package.yml

@@ -42,7 +42,7 @@ jobs:
       # run this locally to update tests durations
       # pytest --cov=lobsterpy --cov-append --splits 1 --group 1 --durations-path ./lobsterpy/TestData/.pytest-split-durations --store-durations
       run: |
-        pytest --cov=lobsterpy --cov-append --splits 5 --group ${{ matrix.split }} --durations-path ./lobsterpy/TestData/.pytest-split-durations
+        pytest --cov=lobsterpy --cov-report term-missing --cov-append  --splits 5 --group ${{ matrix.split }} --durations-path ./lobsterpy/TestData/.pytest-split-durations
     - name: Upload coverage
       uses: actions/upload-artifact@v3
       with:
@@ -68,4 +68,4 @@ jobs:
           continue-on-error: true
           run: |
             coverage combine coverage*/.coverage*
-            coverage report
+            coverage report --show-missing

lobsterpy/cli.py

@@ -13,10 +13,10 @@
 
 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 PlainCohpPlotter, get_style_list, InteractiveCohpPlotter
+from lobsterpy.plotting import PlainCohpPlotter, get_style_list, IcohpDistancePlotter
 
 
 def main() -> None:
@@ -75,6 +75,19 @@ def get_parser() -> argparse.ArgumentParser:
         type=Path,
         help=('path to POTCAR. Default is "POTCAR".'),
     )
+    input_coops_cobis = input_file_group.add_mutually_exclusive_group()
+    input_coops_cobis.add_argument(
+        "--cobis",
+        "--cobi",
+        action="store_true",
+        help="Specifies input file contains COBIS",
+    )
+    input_coops_cobis.add_argument(
+        "--coops",
+        "--coop",
+        action="store_true",
+        help="Specifies input file contains COOPS",
+    )
 
     output_parent = argparse.ArgumentParser(add_help=False)
     output_file_group = output_parent.add_argument_group("Output files")
@@ -271,6 +284,13 @@ def get_parser() -> argparse.ArgumentParser:
         ),
     )
 
+    subparsers.add_parser(
+        "plot-icohps-distances",
+        aliases=["ploticohpsdistances"],
+        parents=[input_parent, plotting_parent],
+        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",
@@ -284,19 +304,6 @@ def get_parser() -> argparse.ArgumentParser:
         type=int,
         help="List of bond numbers, determining COHPs/COBIs/COOPs to include in plot.",
     )
-    plot_coops_cobis = plot_parser.add_mutually_exclusive_group()
-    plot_coops_cobis.add_argument(
-        "--cobis",
-        "--cobi",
-        action="store_true",
-        help="Plot COBIs",
-    )
-    plot_coops_cobis.add_argument(
-        "--coops",
-        "--coop",
-        action="store_true",
-        help="Plot COOPs",
-    )
     plot_grouping = plot_parser.add_mutually_exclusive_group()
     plot_grouping.add_argument(
         "--summed",
@@ -435,6 +442,8 @@ def run(args):
         "automatic-plot-ia",
         "auto-plot-ia",
         "autoplotia",
+        "plot-icohps-distances",
+        "ploticohpsdistances",
     ]:
         style_kwargs = {}
         style_kwargs.update(_user_figsize(args.width, args.height))
@@ -495,7 +504,9 @@ def run(args):
                 filename = filename.with_name(filename.name + ".gz")
             options = {"are_cobis": False, "are_coops": True}
         else:
-            filename = args.cohpcar
+            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(
@@ -673,6 +684,42 @@ def run(args):
                     'please use "--overwrite" if you would like to overwrite existing lobster inputs'
                 )
 
+    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()

lobsterpy/plotting/__init__.py

@@ -12,8 +12,9 @@
 import numpy as np
 from matplotlib import pyplot as plt
 from pkg_resources import resource_filename
+from pymatgen.io.lobster import Icohplist
 from pymatgen.electronic_structure.core import Spin
-from pymatgen.electronic_structure.cohp import Cohp
+from pymatgen.electronic_structure.cohp import Cohp, IcohpCollection
 from pymatgen.electronic_structure.plotter import CohpPlotter
 import plotly.graph_objs as go
 from lobsterpy.cohp.analyze import Analysis
@@ -648,3 +649,113 @@ def _insert_number_of_bonds_in_label(
              bond label with number of bonds inserted
         """
         return label.replace(character, f"{character} {number_of_bonds} x", 1)
+
+
+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):
+        """
+        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
+
+    def add_icohps(self, label, icohpcollection: IcohpCollection):
+        """
+        Adds a 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,
+        }
+
+    def get_plot(
+        self,
+        ax: "matplotlib.axes.Axes | None" = None,
+        style: "matplotlib.plot.style| 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.
+            style: matplotlib style string, if None, will
+                use lobsterpy style by default.
+            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