diff --git a/src/phelel/velph/cli/ph_selfenergy/cmd_ph_selfenergy.py b/src/phelel/velph/cli/ph_selfenergy/cmd_ph_selfenergy.py index 4fdda86..c041bea 100644 --- a/src/phelel/velph/cli/ph_selfenergy/cmd_ph_selfenergy.py +++ b/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 diff --git a/src/phelel/velph/cli/ph_selfenergy/plot/plot_selfenergy.py b/src/phelel/velph/cli/ph_selfenergy/plot/plot_selfenergy.py index f8373c2..2773a93 100644 --- a/src/phelel/velph/cli/ph_selfenergy/plot/plot_selfenergy.py +++ b/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], ".") diff --git a/src/phelel/velph/utils/plot_eigenvalues.py b/src/phelel/velph/utils/plot_eigenvalues.py index df5764b..14020c9 100644 --- a/src/phelel/velph/utils/plot_eigenvalues.py +++ b/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): diff --git a/src/phelel/velph/utils/vasp.py b/src/phelel/velph/utils/vasp.py index 0b7473d..ae697a8 100644 --- a/src/phelel/velph/utils/vasp.py +++ b/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