Merge pull request #13 from phonopy/fix-el_bands

Update velph-el_bands-plot

src/phelel/api_phelel.py

@@ -488,11 +488,19 @@ def _get_phonopy(
         supercell_matrix: Optional[Union[int, float, Sequence, np.ndarray]],
         primitive_matrix: Optional[Union[str, Sequence, np.ndarray]],
     ) -> Phonopy:
+        """Return Phonopy instance.
+
+        Note
+        ----
+        store_dense_svecs=False is necessary to be compatible with code in VASP.
+
+        """
         return Phonopy(
             self._unitcell,
             supercell_matrix=supercell_matrix,
             primitive_matrix=primitive_matrix,
             symprec=self._symprec,
             is_symmetry=self._is_symmetry,
+            store_dense_svecs=False,
             log_level=self._log_level,
         )

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

@@ -51,16 +51,11 @@ def cmd_generate(toml_filename: str):
 def cmd_plot(window: tuple[float, float]):
     """Plot electronic band structure."""
     vaspout_filename_bands = pathlib.Path("el_bands/bands/vaspout.h5")
-    if vaspout_filename_bands.exists():
-        vaspout_filename_dos = pathlib.Path("el_bands/dos/vaspout.h5")
-        if vaspout_filename_bands.exists():
-            click.echo(f'Found "{vaspout_filename_dos}". DOS will be plotted.')
-            plot_el_bandstructures(
-                window,
-                vaspout_filename_bands,
-                vaspout_filename_dos=vaspout_filename_dos,
-            )
-        else:
-            plot_el_bandstructures(window, vaspout_filename_bands)
-    else:
+    vaspout_filename_dos = pathlib.Path("el_bands/dos/vaspout.h5")
+    if not vaspout_filename_bands.exists():
         click.echo(f'"{vaspout_filename_bands}" not found.')
+    if not vaspout_filename_dos.exists():
+        click.echo(f'"{vaspout_filename_dos}" not found.')
+
+    if vaspout_filename_bands.exists() and vaspout_filename_dos.exists():
+        plot_el_bandstructures(window, vaspout_filename_bands, vaspout_filename_dos)

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

@@ -3,7 +3,6 @@
 from __future__ import annotations
 
 import pathlib
-from typing import Optional
 
 import click
 import h5py
@@ -19,7 +18,7 @@
 def plot_el_bandstructures(
     window: tuple[float, float],
     vaspout_filename_bands: pathlib.Path,
-    vaspout_filename_dos: Optional[pathlib.Path] = None,
+    vaspout_filename_dos: pathlib.Path,
     plot_filename: pathlib.Path = pathlib.Path("el_bands/el_bands.pdf"),
 ):
     """Plot electronic band structure.
@@ -37,15 +36,35 @@ 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"][()]
+    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)
+
+    ax0.plot(distances, eigvals[0, :, :], ".k")
+    ax0.hlines(efermi, distances[0], distances[-1], "r")
+    for x in points[1:-1]:
+        ax0.vlines(x, efermi + emin, efermi + emax, "k")
+    ax0.set_xlim(distances[0], distances[-1])
+    ax0.set_xticks(points)
+    ax0.set_xticklabels(labels_at_points)
+    ax0.set_ylim(efermi + emin, efermi + emax)
+    ax0.set_ylabel("E[eV]", fontsize=14)
+    ymin, ymax = ax0.get_ylim()
+
+    print(ymin, ymax)
+    dos, energies, xmax = _get_dos_data(f_h5py_dos, ymin, ymax)
+
+    ax1.plot(dos, energies, "-k")
+    ax1.hlines(efermi, 0, xmax, "r")
+    ax1.set_xlim(0, xmax)
+    ax1.set_ylim(ymin, ymax)
+    ax1.yaxis.tick_right()
 
-    if vaspout_filename_dos:
-        f_h5py_dos = h5py.File(vaspout_filename_dos)
-        _, axs = plt.subplots(1, 2, gridspec_kw={"width_ratios": [3, 1]})
-        _plot_bands(axs[0], f_h5py_bands, window)
-        _plot_dos(axs[1], f_h5py_dos, axs[0])
-    else:
-        _, ax = plt.subplots()
-        _plot_bands(ax, f_h5py_bands, window)
     plt.rcParams["pdf.fonttype"] = 42
     plt.tight_layout()
     plt.savefig(plot_filename)
@@ -54,11 +73,7 @@ def plot_el_bandstructures(
     click.echo(f'Electronic band structure plot was saved in "{plot_filename}".')
 
 
-def _plot_bands(ax, f_h5py, window: tuple[float, float]):
-    emin = window[0]
-    emax = window[1]
-
-    efermi = f_h5py["results"]["electron_dos"]["efermi"][()]
+def _get_bands_data(f_h5py: h5py.File):
     eigvals = f_h5py["results"]["electron_eigenvalues_kpoints_opt"]["eigenvalues"]
 
     reclat = get_reclat_from_vaspout(f_h5py)
@@ -81,23 +96,12 @@ def _plot_bands(ax, f_h5py, window: tuple[float, float]):
         labels, distances, n_segments, nk_per_seg, nk_total
     )
 
-    ax.plot(distances, eigvals[0, :, :], ".k")
-    ax.hlines(efermi, distances[0], distances[-1], "r")
-    for x in points[1:-1]:
-        ax.vlines(x, efermi + emin, efermi + emax, "k")
-    ax.set_xlim(distances[0], distances[-1])
-    ax.set_xticks(points)
-    ax.set_xticklabels(labels_at_points)
-    ax.set_ylim(efermi + emin, efermi + emax)
-    ax.set_ylabel("E[eV]", fontsize=14)
+    return distances, eigvals, points, labels_at_points
 
 
-def _plot_dos(ax, f_h5py, ax_bands):
-    efermi = f_h5py["results"]["electron_dos_kpoints_opt"]["efermi"][()]
+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"][:]
-    ax.plot(dos, energies, "-k")
-    ymin, ymax = ax_bands.get_ylim()
     i_min = 0
     i_max = len(energies)
     for i, val in enumerate(energies):
@@ -107,7 +111,4 @@ def _plot_dos(ax, f_h5py, ax_bands):
             i_max = i
             break
     xmax = dos[i_min : i_max + 1].max() * 1.1
-    ax.hlines(efermi, 0, xmax, "r")
-    ax.set_xlim(0, xmax)
-    ax.set_ylim(ymin, ymax)
-    ax.yaxis.tick_right()
+    return dos, energies, xmax

src/phelel/velph/cli/transport/plot/__init__.py

@@ -72,7 +72,7 @@ def _get_f_h5py_and_plot_filename(
     property_name: str,
     vaspout_filename: pathlib.Path = pathlib.Path("transport/vaspout.h5"),
     plot_filename: Optional[pathlib.Path] = None,
-) -> tuple[h5py._hl.files.File, pathlib.Path]:
+) -> tuple[h5py.File, pathlib.Path]:
     if not vaspout_filename.exists():
         click.echo(f'"{vaspout_filename}" (default path) not found.')
         click.echo("Please specify vaspout.h5 file path.")

src/phelel/velph/cli/transport/plot/plot_selfenergy.py

@@ -7,9 +7,7 @@
 import matplotlib.pyplot as plt
 
 
-def plot_selfenergy(
-    f_h5py: h5py._hl.files.File, plot_filename: str, save_plot: bool = False
-):
+def plot_selfenergy(f_h5py: h5py.File, plot_filename: str, save_plot: bool = False):
     """Plot imaginary part of self-energies.
 
     Number of "self_energy_*" is
@@ -60,7 +58,7 @@ def _plot(ax, selfen):
             )
 
 
-def _show(selfen: h5py._hl.group.Group, index: int):
+def _show(selfen: h5py.Group, index: int):
     """Show self-energy properties.
 
     ['band_start', 'band_stop', 'bks_idx', 'carrier_per_cell',

src/phelel/velph/cli/transport/plot/plot_transport.py

@@ -9,9 +9,7 @@
 import numpy as np
 
 
-def plot_transport(
-    f_h5py: h5py._hl.files.File, plot_filename: str, save_plot: bool = False
-):
+def plot_transport(f_h5py: h5py.File, plot_filename: str, save_plot: bool = False):
     """Plot transport properties.
 
     Number of "transport_*" is
@@ -118,7 +116,7 @@ def _plot(axs: np.ndarray, transports_temps: dict, property_names: tuple):
         ax.set_yticklabels([])
 
 
-def _show(transport: h5py._hl.group.Group, index: int):
+def _show(transport: h5py.Group, index: int):
     """Show transport properties.
 
     ['cell_volume', 'dfermi_tol', 'e_conductivity', 'e_t_conductivity', 'emax',

src/phelel/velph/cli/utils.py

@@ -243,7 +243,7 @@ def choose_cell_in_dict(
     return cell
 
 
-def get_reclat_from_vaspout(fp_vaspout: h5py._hl.files.File):
+def get_reclat_from_vaspout(fp_vaspout: h5py.File):
     """Return reciprocal space basis vectors.
 
     Returns
@@ -255,7 +255,7 @@ def get_reclat_from_vaspout(fp_vaspout: h5py._hl.files.File):
     """
     # Basis vectors in direct space in column vectors
     lattice = np.transpose(
-        fp_vaspout["input"]["poscar"]["lattice_vectors"]
+        fp_vaspout["input"]["poscar"]["lattice_vectors"][:]
         * fp_vaspout["input"]["poscar"]["scale"][()]
     )
     # Basis vectors in reciprocal space in row vectors

test/test_api_phelel.py

@@ -42,11 +42,27 @@ def test_read_phelel_params_hdf5(phelel_CdAs2_111: Phelel):
 
 
 def _compare(filename: pathlib.Path, phe: Phelel):
+    """Assert results.
+
+    shortest_vectors and shortest_vector_multiplicities are included at later
+    versions of Phelel. So, these are not included in old reference files.
+
+    """
     with h5py.File(filename, "r") as f:
         dVdu_ref = f["dVdu"][:]
         dDijdu_ref = f["dDijdu"][:]
 
-    dVdu = cmplx2real(phe.dVdu.dVdu)
-    dDijdu = cmplx2real(phe.dDijdu.dDijdu)
-    np.testing.assert_allclose(dVdu, dVdu_ref, rtol=1e-5, atol=1e-5)
-    np.testing.assert_allclose(dDijdu, dDijdu_ref, rtol=1e-5, atol=1e-5)
+        dVdu = cmplx2real(phe.dVdu.dVdu)
+        dDijdu = cmplx2real(phe.dDijdu.dDijdu)
+        np.testing.assert_allclose(dVdu, dVdu_ref, rtol=1e-5, atol=1e-5)
+        np.testing.assert_allclose(dDijdu, dDijdu_ref, rtol=1e-5, atol=1e-5)
+
+        if "shortest_vectors" in f:
+            shortest_vectors_ref = f["shortest_vectors"][:]
+            multiplicities_ref = f["shortest_vector_multiplicities"][:]
+
+            shortest_vectors, multiplicities = phe.primitive.get_smallest_vectors()
+            np.testing.assert_array_equal(
+                shortest_vectors.shape, shortest_vectors_ref.shape
+            )
+            np.testing.assert_array_equal(multiplicities, multiplicities_ref)