Normalized DOS (#2808)

* complete_dos_normalized

* allow normalized dos

* test

* test

* fixed typo

* also scale pdos

* also scale pdos

* OK don't scale pdos

* tweak doc strings and types

* update plotter, while still allowing scientific malpractice.

Also fixed some linting

* add missing default value

Co-authored-by: Janosh Riebesell <janosh.riebesell@gmail.com>

pymatgen/electronic_structure/dos.py

@@ -192,16 +192,23 @@ class Dos(MSONable):
         Fermi level
     """
 
-    def __init__(self, efermi: float, energies: ArrayLike, densities: Mapping[Spin, ArrayLike]):
+    def __init__(
+        self, efermi: float, energies: ArrayLike, densities: Mapping[Spin, ArrayLike], norm_vol: float | None = None
+    ) -> None:
         """
         Args:
             efermi: Fermi level energy
             energies: A sequences of energies
             densities (dict[Spin: np.array]): representing the density of states for each Spin.
+            norm_vol: The volume used to normalize the densities. Defaults to 1 if None which will not perform any
+                normalization. If not None, the resulting density will have units of states/eV/Angstrom^3, otherwise
+                the density will be in states/eV.
         """
         self.efermi = efermi
         self.energies = np.array(energies)
-        self.densities = {k: np.array(d) for k, d in densities.items()}
+        self.norm_vol = norm_vol
+        vol = norm_vol or 1
+        self.densities = {k: np.array(d) / vol for k, d in densities.items()}
 
     def get_densities(self, spin: Spin | None = None):
         """
@@ -648,17 +655,23 @@ def __init__(
         structure: Structure,
         total_dos: Dos,
         pdoss: Mapping[PeriodicSite, Mapping[Orbital, Mapping[Spin, ArrayLike]]],
-    ):
+        normalize: bool = False,
+    ) -> None:
         """
         Args:
             structure: Structure associated with this particular DOS.
             total_dos: total Dos for structure
             pdoss: The pdoss are supplied as an {Site: {Orbital: {Spin:Densities}}}
+            normalize: Whether to normalize the densities by the volume of the structure.
+                If True, the units of the densities are states/eV/Angstrom^3. Otherwise,
+                the units are states/eV.
         """
+        vol = structure.volume if normalize else None
         super().__init__(
             total_dos.efermi,
             energies=total_dos.energies,
             densities={k: np.array(d) for k, d in total_dos.densities.items()},
+            norm_vol=vol,
         )
         self.pdos = pdoss
         self.structure = structure

pymatgen/electronic_structure/plotter.py

@@ -75,6 +75,7 @@ def __init__(self, zero_at_efermi: bool = True, stack: bool = False, sigma: floa
         self.zero_at_efermi = zero_at_efermi
         self.stack = stack
         self.sigma = sigma
+        self._norm_val = True
         self._doses: dict[
             str, dict[Literal["energies", "densities", "efermi"], float | ArrayLike | dict[Spin, ArrayLike]]
         ] = {}
@@ -84,11 +85,11 @@ def add_dos(self, label: str, dos: Dos) -> None:
         Adds a dos for plotting.
 
         Args:
-            label:
-                label for the DOS. Must be unique.
-            dos:
-                Dos object
+            label: label for the DOS. Must be unique.
+            dos: Dos object
         """
+        if dos.norm_vol is None:
+            self._norm_val = False
         energies = dos.energies - dos.efermi if self.zero_at_efermi else dos.energies
         densities = dos.get_smeared_densities(self.sigma) if self.sigma else dos.densities
         efermi = dos.efermi
@@ -215,7 +216,11 @@ def get_plot(self, xlim=None, ylim=None):
             plt.plot([0, 0], ylim, "k--", linewidth=2)
 
         plt.xlabel("Energies (eV)")
-        plt.ylabel("Density of states")
+
+        if self._norm_val:
+            plt.ylabel("Density of states (states/eV/ų)")
+        else:
+            plt.ylabel("Density of states (states/eV)")
 
         plt.axhline(y=0, color="k", linestyle="--", linewidth=2)
         plt.legend()
@@ -315,7 +320,7 @@ def add_bs(self, bs: BandStructureSymmLine | list[BandStructureSymmLine]) -> Non
 
     def _maketicks(self, plt):
         """
-        utility private method to add ticks to a band structure
+        Utility private method to add ticks to a band structure
         """
         ticks = self.get_ticks()
         # Sanitize only plot the uniq values
@@ -595,6 +600,7 @@ def get_plot(
             smooth (bool or list(bools)): interpolates the bands by a spline cubic.
                 A single bool values means to interpolate all the bandstructure objs.
                 A list of bools allows to select the bandstructure obs to interpolate.
+            vbm_cbm_marker (bool): if True, a marker is added to the vbm and cbm.
             smooth_tol (float) : tolerance for fitting spline to band data.
                 Default is None such that no tolerance will be used.
             smooth_k (int): degree of splines 1<k<5
@@ -752,6 +758,8 @@ def save_plot(self, filename, img_format="eps", ylim=None, zero_to_efermi=True,
             filename: Filename to write to.
             img_format: Image format to use. Defaults to EPS.
             ylim: Specifies the y-axis limits.
+            zero_to_efermi: Automatically the Fermi level as the origin.
+            smooth: Cubic spline interpolation of the bands.
         """
         plt = self.get_plot(ylim=ylim, zero_to_efermi=zero_to_efermi, smooth=smooth)
         plt.savefig(filename, format=img_format)
@@ -840,13 +848,14 @@ def get_ticks_old(self):
 
     def plot_compare(self, other_plotter, legend=True):
         """
-        plot two band structure for comparison. One is in red the other in blue
+        Plot two band structure for comparison. One is in red the other in blue
         (no difference in spins). The two band structures need to be defined
         on the same symmetry lines! and the distance between symmetry lines is
         the one of the band structure used to build the BSPlotter
 
         Args:
-            another band structure object defined along the same symmetry lines
+            other_plotter: Another band structure object defined along the same symmetry lines
+            legend: True to add a legend to the plot
 
         Returns:
             a matplotlib object with both band structures
@@ -887,9 +896,7 @@ def plot_compare(self, other_plotter, legend=True):
         return plt
 
     def plot_brillouin(self):
-        """
-        plot the Brillouin zone
-        """
+        """Plot the Brillouin zone"""
         # get labels and lines
         labels = {}
         for k in self._bs[0].kpoints:
@@ -1142,7 +1149,7 @@ def get_elt_projected_plots(self, zero_to_efermi=True, ylim=None, vbm_cbm_marker
 
     def get_elt_projected_plots_color(self, zero_to_efermi=True, elt_ordered=None):
         """
-        returns a pylab plot object with one plot where the band structure
+        Returns a pylab plot object with one plot where the band structure
         line color depends on the character of the band (along different
         elements). Each element is associated with red, green or blue
         and the corresponding rgb color depending on the character of the band
@@ -2131,7 +2138,7 @@ def orbital_label(list_orbitals):
 
     def _maketicks_selected(self, plt, branches):
         """
-        utility private method to add ticks to a band structure with selected branches
+        Utility private method to add ticks to a band structure with selected branches
         """
         ticks = self.get_ticks()
         distance = []
@@ -2585,12 +2592,14 @@ def _rgbline(ax, k, e, red, green, blue, alpha=1, linestyles="solid"):
     def _get_colordata(bs, elements, bs_projection):
         """
         Get color data, including projected band structures
+
         Args:
             bs: Bandstructure object
             elements: elements (in desired order) for setting to blue, red, green
             bs_projection: None for no projection, "elements" for element projection
 
         Returns:
+            Dictionary representation of color data.
         """
         contribs = {}
         if bs_projection and bs_projection.lower() == "elements":
@@ -2805,6 +2814,7 @@ def plot_seebeck_eff_mass_mu(self, temps=(300,), output="average", Lambda=0.5):
             temps:  list of temperatures of calculated seebeck.
             Lambda: fitting parameter used to model the scattering (0.5 means
                 constant relaxation time).
+
         Returns:
             a matplotlib object
         """
@@ -2865,6 +2875,7 @@ def plot_complexity_factor_mu(self, temps=(300,), output="average", Lambda=0.5):
             temps:  list of temperatures of calculated seebeck and conductivity.
             Lambda: fitting parameter used to model the scattering (0.5 means constant
                     relaxation time).
+
         Returns:
             a matplotlib object
         """
@@ -2917,6 +2928,7 @@ def plot_seebeck_mu(self, temp=600, output="eig", xlim=None):
                 the temperature
             xlim:
                 a list of min and max fermi energy by default (0, and band gap)
+
         Returns:
             a matplotlib object
         """
@@ -3046,6 +3058,7 @@ def plot_seebeck_temp(self, doping="all", output="average"):
                      Specify a list of doping levels if you want to plot only some.
             output: with 'average' you get an average of the three directions
                     with 'eigs' you get all the three directions.
+
         Returns:
             a matplotlib object
         """
@@ -3576,8 +3589,7 @@ def plot_eff_mass_dop(self, temps="all", output="average"):
         return plt
 
     def plot_dos(self, sigma=0.05):
-        """
-        plot dos
+        """Plot dos
 
         Args:
             sigma: a smearing
@@ -3935,6 +3947,7 @@ def plot_fermi_surface(
             If False a non interactive figure will be shown, but it is possible
             to plot other surfaces on the same figure. To make it interactive,
             run mlab.show().
+
     Returns:
         ((mayavi.mlab.figure, mayavi.mlab)): The mlab plotter and an interactive
             figure to control the plot.
@@ -4160,7 +4173,6 @@ def plot_path(line, lattice=None, coords_are_cartesian=False, ax=None, **kwargs)
     Returns:
         matplotlib figure and matplotlib ax
     """
-
     ax, fig, plt = get_ax3d_fig_plt(ax)
 
     if "color" not in kwargs:
@@ -4234,7 +4246,6 @@ def fold_point(p, lattice, coords_are_cartesian=False):
     Returns:
         The Cartesian coordinates folded inside the first Brillouin zone
     """
-
     if coords_are_cartesian:
         p = lattice.get_fractional_coords(p)
     else:
@@ -4354,7 +4365,6 @@ def plot_brillouin_zone(
     Returns:
         matplotlib figure
     """
-
     fig, ax = plot_lattice_vectors(bz_lattice, ax=ax)
     plot_wigner_seitz(bz_lattice, ax=ax)
     if lines is not None:
@@ -4416,13 +4426,14 @@ def plot_ellipsoid(
         kwargs: kwargs passed to the matplotlib function 'plot_wireframe'.
                 Color defaults to blue, rstride and cstride
                 default to 4, alpha defaults to 0.2.
+
     Returns:
         matplotlib figure and matplotlib ax
+
     Example of use:
         fig,ax=plot_wigner_seitz(struct.reciprocal_lattice)
         plot_ellipsoid(hessian,[0.0,0.0,0.0], struct.reciprocal_lattice,ax=ax)
     """
-
     if (not coords_are_cartesian) and lattice is None:
         raise ValueError("coords_are_cartesian False or fold True require the lattice")
 

pymatgen/io/vasp/outputs.py

@@ -692,6 +692,17 @@ def complete_dos(self):
         pdoss = {final_struct[i]: pdos for i, pdos in enumerate(self.pdos)}
         return CompleteDos(self.final_structure, self.tdos, pdoss)
 
+    @property
+    def complete_dos_normalized(self) -> CompleteDos:
+        """
+        A CompleteDos object which incorporates the total DOS and all
+        projected DOS. Normalized by the volume of the unit cell with
+        units of states/eV/unit cell volume.
+        """
+        final_struct = self.final_structure
+        pdoss = {final_struct[i]: pdos for i, pdos in enumerate(self.pdos)}
+        return CompleteDos(self.final_structure, self.tdos, pdoss, normalize=True)
+
     @property
     def hubbards(self):
         """

pymatgen/io/vasp/tests/test_outputs.py

@@ -198,6 +198,14 @@ def test_standard(self):
         self.assertAlmostEqual(pdos0[Orbital.pz][Spin.down][16], 0.0012)
         self.assertEqual(pdos0[Orbital.s][Spin.up].shape, (301,))
 
+        pdos0_norm = vasprun.complete_dos_normalized.pdos[vasprun.final_structure[0]]
+        self.assertAlmostEqual(pdos0_norm[Orbital.s][Spin.up][16], 0.0026)  # the site data should not change
+        self.assertEqual(pdos0_norm[Orbital.s][Spin.up].shape, (301,))
+
+        cdos_norm, cdos = vasprun.complete_dos_normalized, vasprun.complete_dos
+        ratio = np.nanmax(cdos.densities[Spin.up] / cdos_norm.densities[Spin.up])
+        self.assertAlmostEqual(ratio, vasprun.final_structure.volume)  # the site data should not change
+
         filepath2 = self.TEST_FILES_DIR / "lifepo4.xml"
         vasprun_ggau = Vasprun(filepath2, parse_projected_eigen=True, parse_potcar_file=False)
         totalscsteps = sum(len(i["electronic_steps"]) for i in vasprun.ionic_steps)