Added velph ph_selfenergy dump_phono3py

src/phelel/velph/cli/ph_selfenergy/cmd_ph_selfenergy.py

@@ -4,10 +4,16 @@
 
 import click
 import h5py
+import numpy as np
+from numpy.typing import NDArray
 
 from phelel.velph.cli.cmd_root import cmd_root
 from phelel.velph.cli.selfenergy.generate import write_selfenergy_input_files
 from phelel.velph.cli.utils import check_fft
+from phelel.velph.utils.vasp import (
+    convert_ir_kpoints_from_VASP_to_phono3py,
+    read_freqs_and_ph_gammas_from_vaspout_h5,
+)
 
 
 @cmd_root.group("ph_selfenergy")
@@ -54,7 +60,7 @@ def cmd_check_fft(toml_filename: str):
     check_fft(toml_filename, "ph_selfenergy")
 
 
-@cmd_ph_selfenergy.command("dump")
+@cmd_ph_selfenergy.command("dump_phono3py")
 @click.argument(
     "vaspout_filename",
     nargs=1,
@@ -65,13 +71,78 @@ def cmd_check_fft(toml_filename: str):
     "output_filename",
     nargs=1,
     type=click.Path(),
-    default="ph_selfenergy/ph_selfenergy.hdf5",
+    default="ph_selfenergy/phono3py_elph.hdf5",
 )
 @click.help_option("-h", "--help")
-def cmd_dump(vaspout_filename: str, output_filename: str):
-    """Dump ph_selfenergy data to HDF5 file."""
+def cmd_dump_phono3py(vaspout_filename: str, output_filename: str):
+    """Dump ph_selfenergy data to HDF5 file.
+
+    gammas: [ispin, ib, ikpt , nw, temp] -> [ispin, temp, ikpt, ib]
+    freqs: [ikpt, ib] -> [ikpt, ib]
+
+    """
     with h5py.File(vaspout_filename, "r") as f:
-        list(f)
+        id_map, ir_kpoints, weights = _collect_data_from_vaspout(f)
+        freqs_calcs, gammas_calcs, temps_calcs, indices = (
+            read_freqs_and_ph_gammas_from_vaspout_h5(f)
+        )
+        with h5py.File(output_filename, "w") as f_out:
+            for i, freqs, gammas, temps in zip(
+                indices, freqs_calcs, gammas_calcs, temps_calcs, strict=True
+            ):
+                _freqs = freqs[id_map] / (2 * np.pi)  # convert to THz
+                _gammas = np.transpose(gammas[:, :, id_map, 0, :], (0, 3, 2, 1)) / (
+                    2 * np.pi
+                )  # convert to THz
+                f_out.create_dataset(f"frequency_{i}", data=_freqs)
+                f_out.create_dataset(f"gamma_{i}", data=_gammas)
+                f_out.create_dataset(f"temperature_{i}", data=temps)
+                f_out.create_dataset("qpoint", data=ir_kpoints)
+                f_out.create_dataset("weight", data=weights)
+
+            click.echo(f'Dumped ph_selfenergy data to "{output_filename}".')
+
+
+def _collect_data_from_vaspout(f_h5py: h5py.File) -> tuple[NDArray, NDArray, NDArray]:
+    """Collect crystal structure information from vaspout.h5 file."""
+    # [ikpt, 3]
+    ir_kpoints: NDArray = f_h5py[
+        "results/electron_phonon/electrons/eigenvalues/kpoint_coords"
+    ][:]  # type: ignore
+    ir_kpoints_weights: NDArray = f_h5py[
+        "results/electron_phonon/electrons/eigenvalues/kpoints_symmetry_weight"
+    ][:]  # type: ignore
+    lat_scale = f_h5py["results/positions/scale"][()]  # type:ignore
+    # lattice: row vectors
+    lattice: NDArray = f_h5py["results/positions/lattice_vectors"][:] * lat_scale  # type: ignore
+    positions: NDArray = f_h5py["results/positions/position_ions"][:]  # type: ignore
+    number_ion_types: NDArray = f_h5py["results/positions/number_ion_types"][:]  # type: ignore
+    numbers = []
+    for i, nums in enumerate(number_ion_types):
+        numbers += [i + 1] * nums
+    numbers = np.array(numbers, dtype=int)
+    if "basis_vectors" in f_h5py["input/kpoints_elph"]:  # type: ignore
+        # When reading "basis_vectors" in python, the following 3x3 ndarray
+        # is obtained:
+        #   [a_m*_x, a_m*_y, a_m*_z]
+        #   [b_m*_x, b_m*_y, b_m*_z]
+        #   [c_m*_x, c_m*_y, c_m*_z]
+        k_gen_vecs: NDArray = f_h5py["input/kpoints_elph/basis_vectors"][:]  # type:ignore
+    else:
+        k_gen_vecs: NDArray = np.diag(
+            [
+                1.0 / f_h5py[f"input/kpoints_elph/{key}"][()]  # type: ignore
+                for key in ("nkpx", "nkpy", "nkpz")
+            ]
+        )
+    assert f_h5py["input/kpoints/coordinate_space"][()] == b"R"  # type: ignore
+
+    id_map, _, _, _, _ = convert_ir_kpoints_from_VASP_to_phono3py(
+        lattice, positions, numbers, k_gen_vecs, ir_kpoints, ir_kpoints_weights
+    )
+    weights = ir_kpoints_weights[id_map] / np.linalg.det(k_gen_vecs)
+    weights = np.rint(weights).astype(int)
+    return id_map, ir_kpoints[id_map], weights
 
 
 from phelel.velph.cli.ph_selfenergy.plot.cmd_plot import cmd_plot  # noqa: E402, F401

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

@@ -3,9 +3,16 @@
 from __future__ import annotations
 
 import os
+from typing import TYPE_CHECKING
 
 import click
 import h5py
+from numpy.typing import NDArray
+
+from phelel.velph.utils.vasp import read_freqs_and_ph_gammas_from_vaspout_h5
+
+if TYPE_CHECKING:
+    import matplotlib.axes as axes
 
 
 def plot_selfenergy(
@@ -27,32 +34,37 @@ def plot_selfenergy(
     """
     import matplotlib.pyplot as plt
 
-    selfens = {}
-    f_elph = f_h5py["results/electron_phonon/phonons"]
-    for key in f_elph:  # type: ignore
-        if "self_energy_" in key:
-            index = key.split("_")[-1]
-            if index.isdigit():
-                selfens[int(index)] = f_elph[key]  # type: ignore
-
-    assert len(selfens) == int(
-        f_h5py["results/electron_phonon/electrons/self_energy_meta/ncalculators"][()]  # type: ignore
+    freqs_calcs, gammas_calcs, temps_calcs, indices = (
+        read_freqs_and_ph_gammas_from_vaspout_h5(f_h5py)
     )
 
-    if len(selfens) == 1:
-        fig, axs = plt.subplots(1, 1, figsize=(4, 4))
+    n_temp = len(temps_calcs[0])
+    n_spin = gammas_calcs[0].shape[0]
+    n_plots = len(gammas_calcs) * n_spin * n_temp
+    if n_plots == 1:
+        _, axs = plt.subplots(1, 1, figsize=(4, 4))
     else:
-        nrows = len(selfens) // 2
-        fig, axs = plt.subplots(nrows, 2, figsize=(8, 4 * nrows), squeeze=True)
-
-    for i in range(len(selfens)):
-        selfen = selfens[i + 1]
-        if len(selfens) == 1:
-            _axs = axs
-        else:
-            _axs = axs[i]  # type: ignore
-        _plot(_axs, selfen)
-        _axs.set_yscale("log")
+        nrows = n_plots
+        _, axs = plt.subplots(nrows, 1, figsize=(4, 4 * nrows), squeeze=True)
+
+    i_plot = 0
+    for i, (gammas, freqs, temps) in enumerate(
+        zip(gammas_calcs, freqs_calcs, temps_calcs, strict=True)
+    ):
+        for i_spin, gammas_at_spin in enumerate(gammas):
+            for i_temp, temp in enumerate(temps):
+                gammas_at_spin_temp = gammas_at_spin[:, :, :, i_temp]
+                if len(gammas_calcs) * n_spin * n_temp == 1:
+                    ax = axs
+                else:
+                    ax = axs[i_plot]  # type: ignore
+                    i_plot += 1
+                _plot(ax, gammas_at_spin_temp, freqs)
+                ax.set_yscale("log")
+                ax.set_title(
+                    f"$\\tau^{{-1}}$ vs $\\omega$ at T={temp}K spn={i_spin} "
+                    f"({indices[i]}) in 2$\\pi$THz"
+                )
 
     plt.tight_layout()
     if save_plot:
@@ -64,11 +76,8 @@ def plot_selfenergy(
     plt.close()
 
 
-def _plot(ax, selfen):
-    for i_nw in range(selfen["nw"][()]):
-        for selfen_at_T in -selfen["selfen_bubble"][:, i_nw, :, 1].T:
-            ax.plot(
-                selfen["energies"][:, i_nw],
-                selfen_at_T,
-                ".",
-            )
+def _plot(ax: axes.Axes, gammas: NDArray, freqs: NDArray):
+    nw = gammas.shape[2]
+    for i_nw in range(nw):
+        for gammas_at_band, freqs_at_band in zip(gammas, freqs.T, strict=True):
+            ax.plot(freqs_at_band, gammas_at_band[:, i_nw], ".")

src/phelel/velph/utils/plot_eigenvalues.py

@@ -15,7 +15,7 @@
 from scipy.spatial import Voronoi
 
 from phelel.velph.utils.structure import get_symmetry_dataset
-from phelel.velph.utils.vasp import read_crystal_structure_from_h5
+from phelel.velph.utils.vasp import read_crystal_structure_from_vaspout_h5
 
 if TYPE_CHECKING:
     from mpl_toolkits.mplot3d import Axes3D
@@ -145,7 +145,7 @@ def _plot_eigenvalues(
         Chemical potential in eV. If None, the Fermi energy.
 
     """
-    cell = read_crystal_structure_from_h5(f_h5py, "results/positions")
+    cell = read_crystal_structure_from_vaspout_h5(f_h5py, "results/positions")
     sym_dataset = get_symmetry_dataset(cell)
     rotations = [r.T for r in sym_dataset.rotations]
 
@@ -282,8 +282,8 @@ def _get_ax_3D() -> Axes3D:
     fig = plt.figure()
     ax = fig.add_subplot(1, 1, 1, projection="3d")
     ax.grid(False)
-    ax.set_box_aspect([1, 1, 1])
-    return ax
+    ax.set_box_aspect([1, 1, 1])  # type: ignore
+    return ax  # type: ignore
 
 
 def _plot_Brillouin_zone(bz_lattice: NDArray, ax: Axes3D):

src/phelel/velph/utils/vasp.py

@@ -12,6 +12,7 @@
 import spglib
 from numpy.typing import NDArray
 from phono3py.phonon.grid import BZGrid, get_grid_point_from_address, get_ir_grid_points
+from phonopy.physical_units import get_physical_units
 from phonopy.structure.cells import PhonopyAtoms
 from spglib import SpgCell
 
@@ -502,7 +503,9 @@ def read_magmom(magmom: str) -> list[float] | None:
     return VaspIncar().expand(magmom.strip())
 
 
-def read_crystal_structure_from_h5(f_vaspout_h5: h5py.File, group: str) -> PhonopyAtoms:
+def read_crystal_structure_from_vaspout_h5(
+    f_vaspout_h5: h5py.File, group: str
+) -> PhonopyAtoms:
     """Read crystal structure from vaspout.h5."""
     direct = int(f_vaspout_h5[f"{group}/direct_coordinates"][()])  # type: ignore
     scale = f_vaspout_h5[f"{group}/scale"][()]  # type: ignore
@@ -573,7 +576,58 @@ def convert_ir_kpoints_from_VASP_to_phono3py(
     gps = get_grid_point_from_address(ir_addresss, bz_grid.D_diag)
     irgp = ir_grid_map[gps]
     id_map = np.array([np.where(irgp == gp)[0][0] for gp in ir_grid_points], dtype=int)
-    ir_kpoints_weights *= np.linalg.det(mesh)
-    assert (np.abs(ir_grid_weights - ir_kpoints_weights[id_map]) < 1e-8).all()
+    ir_kpoints_weights_vasp = ir_kpoints_weights * np.linalg.det(mesh)
+    assert (np.abs(ir_grid_weights - ir_kpoints_weights_vasp[id_map]) < 1e-8).all()
 
     return id_map, bz_grid, ir_grid_points, ir_grid_weights, ir_grid_map
+
+
+def read_freqs_and_ph_gammas_from_vaspout_h5(
+    f_h5py: h5py.File,
+) -> tuple[list[NDArray], list[NDArray], list[NDArray], list[int]]:
+    """Read phonon frequencies and phonon full-linewidths from vaspout.h5.
+
+    Returns
+    -------
+    freqs_calcs : list[np.ndarray]
+        List of phonon frequencies for each self-energy calculation in 2piTHz.
+        [[ikpt, ib], ...]
+    gammas_calcs : list[np.ndarray]
+        List of imaginary part of phonon self-energies for each self-energy
+        calculation in 2piTHz. Positive sign is chosen.
+        [[ispin, ib, ikpt , nw, temp], ...]
+    indices : list[int]
+        List of self-energy calculation indices. [1, 2, ...] for self_energy_1,
+        self_energy_2, ...
+
+    """
+    THzToEv = get_physical_units().THzToEv
+
+    f_elph = f_h5py["results/electron_phonon/phonons"]
+    indices = []
+    for key in f_elph:  # type: ignore
+        if "self_energy_" in key:
+            # self_energy_1, self_energy_2, ...
+            index = key.split("_")[-1]
+            if index.isdigit():
+                indices.append(int(index))
+
+    gammas_calcs = []
+    freqs_calcs = []
+    temps_calcs = []
+    for index in sorted(indices):
+        selfen = f_elph[f"self_energy_{index}"]  # type: ignore
+        # Imag part of self-energy [ispin, ib, ikpt , nw, temp] in eV
+        selfen_ph: NDArray = selfen["selfen_ph"][:, :, :, :, :, 1]  # type: ignore
+        # [ikpt, ib] in 2piTHz
+        freqs: NDArray = selfen["phonon_freqs_ibz"][:, :]  # type: ignore
+        temps: NDArray = selfen["temps"][:]  # type: ignore
+        gammas_calcs.append(-selfen_ph / (THzToEv / (2 * np.pi)))
+        freqs_calcs.append(freqs)
+        temps_calcs.append(temps)
+
+    assert len(gammas_calcs) == int(
+        f_h5py["results/electron_phonon/electrons/self_energy_meta/ncalculators"][()]  # type: ignore
+    )
+
+    return freqs_calcs, gammas_calcs, temps_calcs, indices