Merge pull request #388 from naik-aakash/update_featurizers

fix sitewise bwdf

Author: naik-aakash

README.md

@@ -92,6 +92,10 @@ If you use any of the following Featurizers, also cite the respective papers:
 * `FeaturizeCharges`: R. Nelson, C. Ertural, P. C. Müller, R. Dronskowski, in Comprehensive Inorganic Chemistry III, *Elsevier*, **2023**, pp. 141–201. [https://doi.org/10.1016/B978-0-12-823144-9.00120-5](https://doi.org/10.1016/B978-0-12-823144-9.00120-5)
 * `FeaturizeIcoxxlist`: V. L. Deringer, W. Zhang, M. Lumeij, S. Maintz, M. Wuttig, R. Mazzarello, R. Dronskowski, *Angewandte Chemie International Edition 2014*, *53*, 10817–10820. [https://doi.org/10.1002/anie.201404223](https://doi.org/10.1002/anie.201404223)
 
+    Cite the following paper as well if using aysmmetry index:
+
+    * F. Belli, E. Zurek, I. Errea, 2025 [https://arxiv.org/abs/2501.14420](https://arxiv.org/abs/2501.14420).
+
 Please cite [pymatgen](https://github.com/materialsproject/pymatgen), [Lobster](https://schmeling.ac.rwth-aachen.de/cohp/index.php?menuID=1), and [ChemEnv](https://doi.org/10.1107/S2052520620007994) correctly as well.
 
 

docs/fundamentals/index.ipynb

@@ -281,8 +281,9 @@
    ]
   },
   {
-   "metadata": {},
    "cell_type": "markdown",
+   "id": "a2f623edafa5efc9",
+   "metadata": {},
    "source": [
     "### ICOXX based features (ICOXX: ICOHP / ICOBI / ICOOP)\n",
     "\n",
@@ -292,9 +293,18 @@
     "\n",
     "In the formula above, $\\delta$ is the Dirac delta function, $r$ is the distance between atoms A and B, and $B_{\\mathrm{AB}}$ is the bond strength between atoms A and B. The bond strength can be ICOHPs, ICOBIs, or ICOOPs. The BWDF is thus a histogram of bond strengths as a function of bond length. One can compute this for entire structure, for each unique atom pair in the structure, per site or per bond label.\n",
     "\n",
-    "Here, from the BWDF mainly statistical features like mean, standard deviation, skewness, kurtosis, weighted mean, and weighted standard deviation are computed.\n"
-   ],
-   "id": "a2f623edafa5efc9"
+    "Here, from the BWDF mainly statistical features like mean, standard deviation, skewness, kurtosis, weighted mean, and weighted standard deviation are computed.\n",
+    "\n",
+    "As proposed by [F. Belli, E. Zurek, I. Errea, 2025](https://arxiv.org/abs/2501.14420), one can also generate the asymmetry index (ASI) of the local bonding environment using ICOBI / ICOHP / ICOOP for a site using following formulation\n",
+    "\n",
+    "$$\n",
+    "\\boldsymbol{V}_{x} = \\frac{1}{B_{x}} \\sum_{\\alpha=1}^{B_{x}} \\mathrm{iCOXX}(x, \\alpha) \\, \\boldsymbol{i}_{x \\alpha}\n",
+    "$$\n",
+    "\n",
+    "$V_x$, is computed for each atom $x$ in the structure.  The vector is the sum of the iCOXX values for all of the interactions between an atom and its  neighboring atoms $\\alpha$ (iCOXX(x, $\\alpha$)) weighted by a unit vector, $i_{xα}$, which denotes the direction of  each interaction. This summation is performed over all neighbors for which the iCOXX values fall  above a user-defined threshold, and divided by the number of interactions considered, $B_x$.\n",
+    "\n",
+    "One can then get statistics of ASI from all sites in the structure."
+   ]
   }
  ],
  "metadata": {

docs/tutorial/commandlineinterface.rst

@@ -228,6 +228,10 @@ command. Below is an example and sample output using this command.
    :file: tutorial_assets/CdF2.html
 
 
+-  ``lobsterpy plot-auto-ia --save-plot-json auto_plot_data.json`` can be used to save underlying plot data from the automatic analysis as a json file. This json file can be then loaded using monty package via ``monty.serialization.loadfn`` function.
+
+Note that ``--save-plot-json`` also works in combination with  ``lobsterpy plot-auto`` subcommand.
+
 2. Plotting of COHPs/COBIs/COOPs
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -274,7 +278,7 @@ You can plot COHPs/COBIs/COOPs from the command line.
 
 .. code:: bash
 
-    lobsterpy plot-bwdf
+    lobsterpy plot-bwdf --plot-negative --sigma 0.15
 
 .. image:: tutorial_assets/example_bwdf.png
 

docs/tutorial/tutorial.ipynb

@@ -1203,7 +1203,7 @@
    "outputs": [],
    "source": [
     "# get the BWDF stats df\n",
-    "batch_icohp.get_df()"
+    "batch_icohp.get_bwdf_df()"
    ]
   },
   {

docs/tutorial/tutorial_assets/example_bwdf.png

@@ -599,6 +599,14 @@ def get_parser() -> argparse.ArgumentParser:
         help="Normalization of the BWDFs. Options: 'formula_units', 'area', 'counts' and 'none'. "
         "Default: 'formula_units'.",
     )
+    bwdf_plotting_group.add_argument(
+        "-plotneg",
+        "--plot-negative",
+        dest="plotneg",
+        action="store_true",
+        default=False,
+        help="If True, will plot -1*ICOHPs. Works only for ICOHPs. Default: True.",
+    )
     bwdf_plotting_group.add_argument(
         "-siteindex",
         "--site-index",
@@ -1092,7 +1100,7 @@ def run(args):
         for bwdf_dict in plot_data:
             pair = next(iter(bwdf_dict.keys()))
             bwdf_plotter.add_bwdf(bwdf=bwdf_dict, label=label)
-            plot = bwdf_plotter.get_plot(sigma=args.sigma, xlim=args.xlim, ylim=args.ylim)
+            plot = bwdf_plotter.get_plot(sigma=args.sigma, xlim=args.xlim, ylim=args.ylim, plot_negative=args.plotneg)
             title = f"{args.title} : {pair}" if args.title else ""
             plot.title(title)
 

src/lobsterpy/cli.py

@@ -850,17 +850,25 @@ class BatchIcoxxlistFeaturizer:
     :param max_length: maximum bond length for BWDF computation
     :param min_length: minimum bond length for BWDF computation
     :param normalization: normalization strategy for BWDF
+    :param n_electrons_scaling: bool indicating if ICOXX values should be scaled by number of electrons
     :param bin_width: bin width for BWDF
     :param bwdf_df_type: Type of BWDF dataframe to generate
 
-        - "binned": Binned BWDF function.
-        - "stats": Statistical features of BWDF function.
+        - "binned": Binned BWDF.
+        - "stats": Statistical features of BWDF.
         - "sorted_bwdf": BWDF values sorted by distances, ascending.
         - "sorted_dists": Distances sorted by BWDF values (either only positive or negative),
           sorted descending by absolute values.
     :param sorted_dists_mode: only applies if bwdf_df_type=="sorted_dists".
-        Corresponds to param "mode" of get_sorted_dist_df, defines whether BWDF values above or
+        Corresponds to the param "mode" of get_sorted_dist_df, defines whether BWDF values above or
         below zero are considered for distance featurization.
+    :param stats_type: type of statistics to compute from BWDF.
+
+        - "atompair": compute stats from unique atom pairs BWDFs.
+        - "site": compute stats from site BWDFs.
+        - "summed": compute stats from structure BWDFs.
+        - "all": concatenated dataframe from `atompair`, `site` and `summed` options.
+
     :read_icobis: bool to state to read ICOBILIST.lobster from the path
     :read_icoops: bool to state to read ICOOPLIST.lobster from the path
     :param n_jobs: number of parallel processes to run
@@ -871,9 +879,11 @@ def __init__(
         self,
         path_to_lobster_calcs: str | Path,
         normalization: Literal["formula_units", "area", "counts", "none"] = "formula_units",
+        n_electrons_scaling: bool = False,
         bin_width: float = 0.02,
         bwdf_df_type: Literal["binned", "stats", "sorted_bwdf", "sorted_dists"] = "stats",
         sorted_dists_mode: Literal["positive", "negative"] = "negative",
+        stats_type: Literal["atompair", "site", "summed", "all"] = "all",
         interactions_tol: float = 1e-3,
         max_length: float = 6.0,
         min_length: float = 0.0,
@@ -888,6 +898,8 @@ def __init__(
         :param max_length: maximum bond length for BWDF computation
         :param min_length: minimum bond length for BWDF computation
         :param normalization: normalization strategy for BWDF
+        :param n_electrons_scaling: bool indicating if ICOXX values should be scaled by number of electrons.
+            Only for testing purposes. Should not affect the results in any meaningful way.
         :param bin_width: bin width for BWDF
         :param bwdf_df_type: Type of BWDF dataframe to generate
 
@@ -899,59 +911,76 @@ def __init__(
         :param sorted_dists_mode: only applies if bwdf_df_type=="sorted_dists".
             Corresponds to param "mode" of get_sorted_dist_df, defines whether BWDF values above or
             below zero are considered for distance featurization.
+        :param stats_type: type of BWDF stats to be computed. Only applies if bwdf_df_type=="stats".
+
+            - "atompair": compute stats from unique atom pairs BWDFs.
+            - "site": compute stats from site BWDFs.
+            - "summed": compute stats from structure BWDFs.
+            - "all": concatenated dataframe from `atompair`, `site` and `summed` options.
         :param interactions_tol: tolerance for interactions
         :param read_icobis: bool to state to read ICOBILIST.lobster from the path
         :param read_icoops: bool to state to read ICOOPLIST.lobster from the path
         :param n_jobs: number of parallel processes to run
         """
         self.path_to_lobster_calcs = path_to_lobster_calcs
         self.normalization = normalization
+        self.n_electrons_scaling = n_electrons_scaling
         self.max_length = max_length
         self.min_length = min_length
         self.bin_width = bin_width
         self.interactions_tol = interactions_tol
         self.bwdf_df_type = bwdf_df_type
         self.sorted_dists_mode = sorted_dists_mode
+        self.stats_type = stats_type
         self.read_icobis = read_icobis
         self.read_icoops = read_icoops
         self.n_jobs = n_jobs
 
-    def _get_icoxxlist_bwdf_df(self, path_to_lobster_calc: str | Path) -> pd.DataFrame:
+    def _get_featurizer_obj(self, path_to_lobster_calc: str | Path) -> FeaturizeIcoxxlist:
         """
-        Featurize ICOXXLIST data using FeaturizeCOXX.
+        Set up the FeaturizeIcoxxlist object based on the provided parameters.
 
         :param path_to_lobster_calc: path to root LOBSTER calculation directory
 
         Returns:
-            A pandas dataframe with computed ICOXXLIST moment features
+            A FeaturizeIcoxxlist object
         """
         if self.read_icobis:
             file_paths = get_file_paths(
                 path_to_lobster_calc=path_to_lobster_calc,
-                requested_files=["structure", "icobilist"],
+                requested_files=["structure", "icobilist", "grosspop"]
+                if self.n_electrons_scaling
+                else ["structure", "icobilist"],
             )
             feat_icoxx = FeaturizeIcoxxlist(
                 path_to_icoxxlist=file_paths.get("icobilist"),
                 path_to_structure=file_paths.get("structure"),
+                path_to_grosspop=file_paths.get("grosspop", None),
                 bin_width=self.bin_width,
                 interactions_tol=self.interactions_tol,
                 normalization=self.normalization,
+                n_electrons_scaling=self.n_electrons_scaling,
                 max_length=self.max_length,
                 min_length=self.min_length,
                 are_cobis=self.read_icobis,
                 are_coops=self.read_icoops,
             )
+
         elif self.read_icoops:
             file_paths = get_file_paths(
                 path_to_lobster_calc=path_to_lobster_calc,
-                requested_files=["structure", "icooplist"],
+                requested_files=["structure", "icooplist", "grosspop"]
+                if self.n_electrons_scaling
+                else ["structure", "icooplist"],
             )
             feat_icoxx = FeaturizeIcoxxlist(
                 path_to_icoxxlist=file_paths.get("icooplist"),
                 path_to_structure=file_paths.get("structure"),
+                path_to_grosspop=file_paths.get("grosspop", None),
                 bin_width=self.bin_width,
                 interactions_tol=self.interactions_tol,
                 normalization=self.normalization,
+                n_electrons_scaling=self.n_electrons_scaling,
                 max_length=self.max_length,
                 min_length=self.min_length,
                 are_cobis=self.read_icobis,
@@ -960,29 +989,59 @@ def _get_icoxxlist_bwdf_df(self, path_to_lobster_calc: str | Path) -> pd.DataFra
         else:
             file_paths = get_file_paths(
                 path_to_lobster_calc=path_to_lobster_calc,
-                requested_files=["structure", "icohplist"],
+                requested_files=["structure", "icohplist", "grosspop"]
+                if self.n_electrons_scaling
+                else ["structure", "icohplist"],
             )
             feat_icoxx = FeaturizeIcoxxlist(
                 path_to_icoxxlist=file_paths.get("icohplist"),
                 path_to_structure=file_paths.get("structure"),
+                path_to_grosspop=file_paths.get("grosspop", None),
                 bin_width=self.bin_width,
                 interactions_tol=self.interactions_tol,
                 normalization=self.normalization,
+                n_electrons_scaling=self.n_electrons_scaling,
                 max_length=self.max_length,
                 min_length=self.min_length,
                 are_cobis=self.read_icobis,
                 are_coops=self.read_icoops,
             )
 
+        return feat_icoxx
+
+    def _get_icoxxlist_bwdf_df(self, path_to_lobster_calc: str | Path) -> pd.DataFrame:
+        """
+        Featurize ICOXXLIST data using  FeaturizeIcoxxlist.
+
+        :param path_to_lobster_calc: path to root LOBSTER calculation directory
+
+        Returns:
+            A pandas dataframe with computed ICOXXLIST moment features
+        """
+        feat_icoxx = self._get_featurizer_obj(path_to_lobster_calc)
+
         if self.bwdf_df_type == "binned":
             return feat_icoxx.get_binned_bwdf_df()
         if self.bwdf_df_type == "sorted_bwdf":
             return feat_icoxx.get_sorted_bwdf_df()
         if self.bwdf_df_type == "sorted_dists":
             return feat_icoxx.get_sorted_dist_df(mode=self.sorted_dists_mode)
-        return feat_icoxx.get_stats_df()
+        return feat_icoxx.get_stats_df(stats_type=self.stats_type)
 
-    def get_df(self) -> pd.DataFrame:
+    def _get_asymmetry_index_df(self, path_to_lobster_calc: str | Path) -> pd.DataFrame:
+        """
+        Generate a pandas dataframe with asymmetry index stats features.
+
+        :param path_to_lobster_calc: path to root LOBSTER calculation directory
+
+        Returns:
+            A pandas dataframe with asymmetry index stats features as columns.
+        """
+        feat_icoxx = self._get_featurizer_obj(path_to_lobster_calc)
+
+        return feat_icoxx.get_asymmetry_index_stats_df()
+
+    def get_bwdf_df(self) -> pd.DataFrame:
         """
         Generate a pandas dataframe with BWDF for all calcs.
 
@@ -1013,3 +1072,32 @@ def get_df(self) -> pd.DataFrame:
         df_icoxxlist.sort_index(inplace=True)  # noqa: PD002
 
         return df_icoxxlist
+
+    def get_asymmetry_index_df(self) -> pd.DataFrame:
+        """
+        Generate a pandas dataframe with site asymmetry index stats for all calcs.
+
+        Returns:
+            A pandas dataframe with site asymmetry index stats as columns.
+        """
+        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))
+        ]
+        row = []
+        with (
+            mp.Pool(processes=self.n_jobs, maxtasksperchild=1) as pool,
+            tqdm(total=len(paths), desc="Generating ASI from ICOXXLIST") as pbar,
+        ):
+            for _, result in enumerate(pool.imap_unordered(self._get_asymmetry_index_df, paths, chunksize=1)):
+                pbar.update()
+                row.append(result)
+
+        df_icoxxlist = pd.concat(row)
+
+        df_icoxxlist.sort_index(inplace=True)  # noqa: PD002
+
+        return df_icoxxlist