Merge pull request #42 from phonopy/el_dos

Support reading non-KPOINTS_OPT DOS results in vaspout.h5

src/phelel/velph/cli/el_bands/plot.py

@@ -6,6 +6,7 @@
 
 import click
 import h5py
+import numpy as np
 
 from phelel.velph.cli.utils import (
     get_distances_along_BZ_path,
@@ -38,13 +39,31 @@ def plot_el_bandstructures(
 
     f_h5py_bands = h5py.File(vaspout_filename_bands)
     f_h5py_dos = h5py.File(vaspout_filename_dos)
-    efermi = f_h5py_dos["results"]["electron_dos_kpoints_opt"]["efermi"][()]
+    if "results" not in f_h5py_bands:
+        raise ValueError(
+            f"No electronic band structure results found in {vaspout_filename_bands}."
+        )
+    if "results" not in f_h5py_dos:
+        raise ValueError(f"No electronic DOS results found in {vaspout_filename_dos}.")
+    if "electron_dos_kpoints_opt" in f_h5py_dos["results"]:
+        f_h5py_dos_results = f_h5py_dos["results/electron_dos_kpoints_opt"]
+    elif "electron_dos" in f_h5py_dos["results"]:
+        f_h5py_dos_results = f_h5py_dos["results/electron_dos"]
+    else:
+        raise ValueError("No electron DOS data found in vaspout.h5.")
+
+    efermi = f_h5py_dos_results["efermi"][()]
     emin = window[0]
     emax = window[1]
     _, axs = plt.subplots(1, 2, gridspec_kw={"width_ratios": [3, 1]})
     ax0, ax1 = axs
 
-    distances, eigvals, points, labels_at_points = _get_bands_data(f_h5py_bands)
+    reclat = get_reclat_from_vaspout(f_h5py_bands)
+    distances, eigvals, points, labels_at_points = _get_bands_data(
+        reclat,
+        f_h5py_bands["results/electron_eigenvalues_kpoints_opt"],
+        f_h5py_bands["input/kpoints_opt"],
+    )
 
     ax0.plot(distances, eigvals[0, :, :], "ok", markersize=1)
     ax0.hlines(efermi, distances[0], distances[-1], "r", linewidth=1)
@@ -57,7 +76,7 @@ def plot_el_bandstructures(
     ax0.set_ylabel("E[eV]", fontsize=14)
     ymin, ymax = ax0.get_ylim()
 
-    dos, energies, xmax = _get_dos_data(f_h5py_dos, ymin, ymax)
+    dos, energies, xmax = _get_dos_data(f_h5py_dos_results, ymin, ymax)
 
     ax1.plot(dos, energies, "-k", linewidth=1)
     ax1.hlines(efermi, 0, xmax, "r", linewidth=1)
@@ -73,20 +92,18 @@ def plot_el_bandstructures(
     click.echo(f'Electronic band structure plot was saved in "{plot_filename}".')
 
 
-def _get_bands_data(f_h5py: h5py.File):
-    eigvals = f_h5py["results"]["electron_eigenvalues_kpoints_opt"]["eigenvalues"]
+def _get_bands_data(
+    reclat: np.ndarray, f_h5py_bands_results: h5py.Group, f_h5py_bands_input: h5py.Group
+):
+    eigvals = f_h5py_bands_results["eigenvalues"]
 
-    reclat = get_reclat_from_vaspout(f_h5py)
     # k-points in reduced coordinates
-    kpoint_coords = f_h5py["results"]["electron_eigenvalues_kpoints_opt"][
-        "kpoint_coords"
-    ]
+    kpoint_coords = f_h5py_bands_results["kpoint_coords"]
     # Special point labels
     labels = [
-        label.decode("utf-8")
-        for label in f_h5py["input"]["kpoints_opt"]["labels_kpoints"][:]
+        label.decode("utf-8") for label in f_h5py_bands_input["labels_kpoints"][:]
     ]
-    nk_per_seg = f_h5py["input"]["kpoints_opt"]["number_kpoints"][()]
+    nk_per_seg = f_h5py_bands_input["number_kpoints"][()]
     nk_total = len(kpoint_coords)
     k_cart = kpoint_coords @ reclat
     n_segments = nk_total // nk_per_seg
@@ -99,9 +116,9 @@ def _get_bands_data(f_h5py: h5py.File):
     return distances, eigvals, points, labels_at_points
 
 
-def _get_dos_data(f_h5py: h5py.File, ymin: float, ymax: float):
-    dos = f_h5py["results"]["electron_dos_kpoints_opt"]["dos"][0, :]
-    energies = f_h5py["results"]["electron_dos_kpoints_opt"]["energies"][:]
+def _get_dos_data(f_h5py_dos_results: h5py.Group, ymin: float, ymax: float):
+    dos = f_h5py_dos_results["dos"][0, :]
+    energies = f_h5py_dos_results["energies"][:]
     i_min = 0
     i_max = len(energies)
     for i, val in enumerate(energies):

test/velph/cli/el_bands/plot/test_plot.py

@@ -7,6 +7,7 @@
 import pytest
 
 from phelel.velph.cli.el_bands.plot import _get_bands_data, _get_dos_data
+from phelel.velph.cli.utils import get_reclat_from_vaspout
 
 cwd = pathlib.Path(__file__).parent
 
@@ -21,9 +22,16 @@ def test_velph_el_bands_plot_TiNiSn():
     assert vaspout_filename_bands.exists()
     f_h5py_bands = h5py.File(vaspout_filename_bands)
     f_h5py_dos = h5py.File(vaspout_filename_dos)
-    distances, eigvals, points, labels_at_points = _get_bands_data(f_h5py_bands)
+    reclat = get_reclat_from_vaspout(f_h5py_bands)
+    distances, eigvals, points, labels_at_points = _get_bands_data(
+        reclat,
+        f_h5py_bands["results/electron_eigenvalues_kpoints_opt"],
+        f_h5py_bands["input/kpoints_opt"],
+    )
     ymin, ymax = 3.575980267703933, 17.575980267703933
-    dos, energies, xmax = _get_dos_data(f_h5py_dos, ymin, ymax)
+    dos, energies, xmax = _get_dos_data(
+        f_h5py_dos["results/electron_dos_kpoints_opt"], ymin, ymax
+    )
 
     assert len(distances) == 306
     assert pytest.approx(distances[100], 1e-5) == 1.421887803385511