From 3904711f351a47b1117543af206e950c56ff9c4c Mon Sep 17 00:00:00 2001 From: Atsushi Togo Date: Mon, 23 Feb 2026 12:11:27 +0900 Subject: [PATCH] Refactor velph's utility functions --- src/phelel/velph/cli/el_bands/plot.py | 48 ++--- src/phelel/velph/cli/init/init.py | 18 +- src/phelel/velph/cli/ph_bands/plot.py | 28 +-- src/phelel/velph/cli/utils.py | 222 +-------------------- src/phelel/velph/utils/plot_eigenvalues.py | 2 +- src/phelel/velph/utils/structure.py | 137 +++++++++++++ src/phelel/velph/utils/vasp.py | 99 +++++++++ test/velph/cli/el_bands/plot/test_plot.py | 14 +- test/velph/cli/init/test_cmd_init.py | 2 +- test/velph/cli/test_utils.py | 3 +- 10 files changed, 299 insertions(+), 274 deletions(-) create mode 100644 src/phelel/velph/utils/structure.py diff --git a/src/phelel/velph/cli/el_bands/plot.py b/src/phelel/velph/cli/el_bands/plot.py index e24c447..c76e374 100644 --- a/src/phelel/velph/cli/el_bands/plot.py +++ b/src/phelel/velph/cli/el_bands/plot.py @@ -7,12 +7,9 @@ import click import h5py import numpy as np +from numpy.typing import NDArray -from phelel.velph.cli.utils import ( - get_distances_along_BZ_path, - get_reclat_from_vaspout, - get_special_points, -) +from phelel.velph.utils.vasp import get_bands_data, get_reclat_from_vaspout def plot_el_bandstructures( @@ -48,14 +45,14 @@ def plot_el_bandstructures( ) 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"]: + if "electron_dos_kpoints_opt" in f_h5py_dos["results"]: # type: ignore f_h5py_dos_results = f_h5py_dos["results/electron_dos_kpoints_opt"] - elif "electron_dos" in f_h5py_dos["results"]: + elif "electron_dos" in f_h5py_dos["results"]: # type: ignore 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"][()] + efermi: float = f_h5py_dos_results["efermi"][()] # type: ignore emin = window[0] emax = window[1] _, axs = plt.subplots(1, 2, gridspec_kw={"width_ratios": [3, 1]}) @@ -64,8 +61,8 @@ def plot_el_bandstructures( 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"], + f_h5py_bands["results/electron_eigenvalues_kpoints_opt"], # type: ignore + f_h5py_bands["input/kpoints_opt"], # type: ignore ) lines = ["-k", "--k"] @@ -81,7 +78,9 @@ 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_results, ymin, ymax) + dos: NDArray = f_h5py_dos_results["dos"][:] # type: ignore + energies: NDArray = f_h5py_dos_results["energies"][:] # type: ignore + xmax = _get_dos_data(dos, energies, ymin, ymax) lines = ["-k", "--k"] for i, dos_spin in enumerate(dos): @@ -104,32 +103,27 @@ def plot_el_bandstructures( def _get_bands_data( - reclat: np.ndarray, f_h5py_bands_results: h5py.Group, f_h5py_bands_input: h5py.Group + reclat: NDArray, f_h5py_bands_results: h5py.Group, f_h5py_bands_input: h5py.Group ): - eigvals = f_h5py_bands_results["eigenvalues"][:] + eigvals: NDArray = f_h5py_bands_results["eigenvalues"][:] # type: ignore # k-points in reduced coordinates - kpoint_coords = f_h5py_bands_results["kpoint_coords"] + kpoint_coords: NDArray = f_h5py_bands_results["kpoint_coords"][:] # type: ignore # Special point labels labels = [ - label.decode("utf-8") for label in f_h5py_bands_input["labels_kpoints"][:] + label.decode("utf-8") + for label in f_h5py_bands_input["labels_kpoints"][:] # type: ignore ] - 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 - assert n_segments * nk_per_seg == nk_total - distances = get_distances_along_BZ_path(nk_total, n_segments, nk_per_seg, k_cart) - points, labels_at_points = get_special_points( - labels, distances, n_segments, nk_per_seg, nk_total + nk_per_seg: int = f_h5py_bands_input["number_kpoints"][()] # type: ignore + + distances, points, labels_at_points = get_bands_data( + kpoint_coords, reclat, nk_per_seg, labels ) return distances, eigvals, points, labels_at_points -def _get_dos_data(f_h5py_dos_results: h5py.Group, ymin: float, ymax: float): - dos = f_h5py_dos_results["dos"][:] - energies = f_h5py_dos_results["energies"][:] +def _get_dos_data(dos: NDArray, energies: NDArray, ymin: float, ymax: float) -> float: i_min = 0 i_max = len(energies) for i, val in enumerate(energies): @@ -154,4 +148,4 @@ def _get_dos_data(f_h5py_dos_results: h5py.Group, ymin: float, ymax: float): dos = np.where(dos > 10000, 0, dos) xmax = dos[:, i_min : i_max + 1].max() * 1.1 - return dos, energies, xmax + return xmax diff --git a/src/phelel/velph/cli/init/init.py b/src/phelel/velph/cli/init/init.py index 2d1d5fe..fc0f928 100644 --- a/src/phelel/velph/cli/init/init.py +++ b/src/phelel/velph/cli/init/init.py @@ -35,12 +35,14 @@ VelphFilePaths, VelphInitOptions, VelphInitParams, +) +from phelel.velph.templates import default_template_dict +from phelel.velph.utils.structure import ( generate_standardized_cells, get_primitive_cell, get_reduced_cell, get_symmetry_dataset, ) -from phelel.velph.templates import default_template_dict from phelel.velph.utils.vasp import CutoffToFFTMesh, VaspIncar from phelel.version import __version__ @@ -549,6 +551,20 @@ def _get_cells( click.echo(f" Reference space-group-type: {spg_type.international_short}") if symmetrize_cell: + if isinstance(sym_dataset, SpglibDataset): + click.echo( + "Crystal structure was standardized based on space-group-type " + f"{sym_dataset.international}." + ) + elif isinstance(sym_dataset, SpglibMagneticDataset): + click.echo( + "Crystal structure was standardized based on magnetic-space-group-type " + f"UNI No.{sym_dataset.uni_number}." + ) + else: + raise ValueError( + "sym_dataset must be SpglibDataset or SpglibMagneticDataset." + ) unitcell, _primitive, tmat = generate_standardized_cells( sym_dataset, tolerance=tolerance ) diff --git a/src/phelel/velph/cli/ph_bands/plot.py b/src/phelel/velph/cli/ph_bands/plot.py index 58352b6..a5372e6 100644 --- a/src/phelel/velph/cli/ph_bands/plot.py +++ b/src/phelel/velph/cli/ph_bands/plot.py @@ -1,16 +1,15 @@ """Implementation of velph-ph_bands-plot.""" +from __future__ import annotations + import pathlib import click import h5py import numpy as np +from numpy.typing import NDArray -from phelel.velph.cli.utils import ( - get_distances_along_BZ_path, - get_reclat_from_vaspout, - get_special_points, -) +from phelel.velph.utils.vasp import get_bands_data, get_reclat_from_vaspout def plot_ph_bandstructures( @@ -37,24 +36,17 @@ def plot_ph_bandstructures( import matplotlib.pyplot as plt f = h5py.File(vaspout_filename) - eigvals = f["results"]["phonons"]["eigenvalues"][:] # phonon eigenvalues + eigvals: NDArray = f["results/phonons/eigenvalues"][:] # type: ignore if use_ordinary_frequency: eigvals /= 2 * np.pi omega_max = 1.1 * eigvals.max() reclat = get_reclat_from_vaspout(f) - labels = [ - label.decode("utf-8") for label in f["input"]["qpoints"]["labels_kpoints"][:] - ] - nk_per_seg = f["input"]["qpoints"]["number_kpoints"][()] - kpoint_coords = f["results"]["phonons"]["kpoint_coords"] - nk_total = len(kpoint_coords) - k_cart = kpoint_coords @ reclat - n_segments = nk_total // nk_per_seg - assert n_segments * nk_per_seg == nk_total - distances = get_distances_along_BZ_path(nk_total, n_segments, nk_per_seg, k_cart) - points, labels_at_points = get_special_points( - labels, distances, n_segments, nk_per_seg, nk_total + labels = [label.decode("utf-8") for label in f["input/qpoints/labels_kpoints"][:]] # type: ignore + nk_per_seg: int = f["input/qpoints/number_kpoints"][()] # type: ignore + kpoint_coords: NDArray = f["results/phonons/kpoint_coords"][:] # type: ignore + distances, points, labels_at_points = get_bands_data( + kpoint_coords, reclat, nk_per_seg, labels ) _, ax = plt.subplots() diff --git a/src/phelel/velph/cli/utils.py b/src/phelel/velph/cli/utils.py index af11f6c..e6eeb7b 100644 --- a/src/phelel/velph/cli/utils.py +++ b/src/phelel/velph/cli/utils.py @@ -12,19 +12,12 @@ from typing import Any, Iterator, Literal import click -import h5py import numpy as np -import spglib from numpy.typing import NDArray from phono3py.phonon.grid import BZGrid -from phonopy.interface.vasp import VasprunxmlExpat, sort_positions_by_symbols +from phonopy.interface.vasp import VasprunxmlExpat from phonopy.physical_units import get_physical_units -from phonopy.structure.atoms import PhonopyAtoms, get_atomic_data, parse_cell_dict -from phonopy.structure.cells import ( - get_primitive, - get_primitive_matrix_by_centring, - get_reduced_bases, -) +from phonopy.structure.atoms import PhonopyAtoms, parse_cell_dict from phonopy.structure.symmetry import symmetrize_borns_and_epsilon from spglib import SpglibDataset, SpglibMagneticDataset @@ -33,6 +26,7 @@ get_sge_scheduler_script, get_slurm_scheduler_script, ) +from phelel.velph.utils.structure import get_symmetry_dataset from phelel.velph.utils.vasp import VaspIncar, VaspKpoints @@ -338,216 +332,6 @@ def choose_cell_in_dict( return cell -def get_reclat_from_vaspout(fp_vaspout: h5py.File): - """Return reciprocal space basis vectors. - - Returns - ------- - reclat : np.ndarray - Reciprocal basis vectors in row vectors. - shape=(3, 3) - - """ - # Basis vectors in direct space in column vectors - lattice = np.transpose( - fp_vaspout["input"]["poscar"]["lattice_vectors"][:] # type: ignore - * fp_vaspout["input"]["poscar"]["scale"][()] # type: ignore - ) - # Basis vectors in reciprocal space in row vectors - reclat = 2 * np.pi * np.linalg.inv(lattice) - return reclat - - -def get_distances_along_BZ_path(nk_total, n_segments, nk_per_seg, k_cart): - """Measure distances of points from origin along paths. - - Returns - ------- - distances : np.ndarray - Distances of points from origin along BZ-paths. - shape=(nk_total,) - - """ - distances = np.zeros(nk_total) - count = 0 - for _ in range(n_segments): - for i_pts in range(nk_per_seg): - # Treatment of jump between equivalent points on BZ boundary - if i_pts == 0: - delta_dist = 0 - else: - delta_dist = np.linalg.norm(k_cart[count, :] - k_cart[count - 1, :]) - distances[count] = distances[count - 1] + delta_dist - count += 1 - return distances - - -def get_special_points( - labels: list[str], - distances: list[float], - n_segments: int, - nk_per_seg: int, - nk_total: int, -) -> tuple[list, list]: - """Plot special points at vertical lines and labels.""" - # Left most - points = [] - labels_at_points = [] - - labels_at_points.append(labels[0][0]) - points.append(distances[0]) - - count = 0 - for i_seg in range(n_segments): - for _ in range(nk_per_seg): - count += 1 - if count != nk_total: - points.append(distances[count]) - if labels[i_seg * 2 + 1] == labels[i_seg * 2 + 2]: - labels_at_points.append(labels[i_seg * 2 + 1][0]) - else: - labels_at_points.append( - f"{labels[i_seg * 2 + 1][0]}|{labels[i_seg * 2 + 2][0]}" - ) - - labels_at_points.append(labels[-1][0]) - points.append(distances[-1]) - - return points, labels_at_points - - -def get_symmetry_dataset( - cell: PhonopyAtoms, tolerance: float = 1e-5 -) -> SpglibDataset | SpglibMagneticDataset: - """Return spglib symmetry dataset.""" - if cell.magnetic_moments is None: - dataset = spglib.get_symmetry_dataset(cell.totuple(), symprec=tolerance) - else: - dataset = spglib.get_magnetic_symmetry_dataset( - cell.totuple(), symprec=tolerance - ) - - assert dataset is not None - - return dataset - - -def get_primitive_cell( - cell: PhonopyAtoms, - sym_dataset: SpglibDataset | SpglibMagneticDataset, - tolerance: float = 1e-5, -) -> tuple[PhonopyAtoms, NDArray]: - """Return primitive cell and transformation matrix. - - This primitive cell is generated from the input cell without - rigid rotation in contrast to `_get_standardized_unitcell`. - - """ - tmat = sym_dataset.transformation_matrix - pmat = get_primitive_matrix_from_dataset(sym_dataset) - total_tmat = np.array(np.linalg.inv(tmat) @ pmat, dtype="double", order="C") - - return ( - get_primitive(cell, primitive_matrix=total_tmat, symprec=tolerance), - total_tmat, - ) - - -def get_primitive_matrix_from_dataset( - sym_dataset: SpglibDataset | SpglibMagneticDataset, -) -> NDArray: - """Return primitive matrix from symmetry dataset.""" - if isinstance(sym_dataset, SpglibDataset): - centring = sym_dataset.international[0] - else: - spg_type = spglib.get_spacegroup_type(sym_dataset.hall_number) - assert spg_type is not None - centring = spg_type.international_short[0] - return get_primitive_matrix_by_centring(centring) - - -def get_reduced_cell( - cell: PhonopyAtoms, method="niggli", tolerance: float = 1e-5 -) -> PhonopyAtoms: - """Return a reduced cell of input cell.""" - reduced_lattice = get_reduced_bases(cell.cell, method=method, tolerance=tolerance) - assert reduced_lattice is not None - reduced_positions = cell.scaled_positions @ ( - cell.cell @ np.linalg.inv(reduced_lattice) - ) - reduced_positions = reduced_positions - np.rint(reduced_positions) - reduced_positions[:, :] += np.where(reduced_positions < 0, 1, 0) - reduced_cell = cell.copy() - reduced_cell.cell = reduced_lattice - reduced_cell.scaled_positions = reduced_positions - return reduced_cell - - -def generate_standardized_cells( - sym_dataset: SpglibDataset | SpglibMagneticDataset, - tolerance: float = 1e-5, -) -> tuple[PhonopyAtoms, PhonopyAtoms, NDArray]: - """Return standardized unit cell and primitive cell.""" - if isinstance(sym_dataset, SpglibDataset): - click.echo( - "Crystal structure was standardized based on space-group-type " - f"{sym_dataset.international}." - ) - else: - click.echo( - "Crystal structure was standardized based on magnetic-space-group-type " - f"UNI No.{sym_dataset.uni_number}." - ) - convcell = get_standardized_unitcell(sym_dataset) - pmat = get_primitive_matrix_from_dataset(sym_dataset) - if (np.abs(pmat - np.eye(3)) < 1e-8).all(): - primitive = convcell - else: - primitive = get_primitive(convcell, primitive_matrix=pmat, symprec=tolerance) - - return convcell, primitive, pmat - - -def get_standardized_unitcell( - dataset: SpglibDataset | SpglibMagneticDataset, -) -> PhonopyAtoms: - """Return conventional unit cell. - - This conventional unit cell can include rigid rotation with respect to - input unit cell for which symmetry was analized. - - Parameters - ---------- - cell : PhonopyAtoms - Input cell. - dataset : SpgliDataset - Symmetry dataset of spglib. - - Returns - ------- - PhonopyAtoms - Convetional unit cell. - - """ - std_positions = dataset.std_positions - std_types = dataset.std_types - _, _, _, perm = sort_positions_by_symbols(std_types, std_positions) - atom_data = get_atomic_data().atom_data - if isinstance(dataset, SpglibDataset): - return PhonopyAtoms( - cell=dataset.std_lattice, - scaled_positions=std_positions[perm], - symbols=[atom_data[n][1] for n in std_types[perm]], - ) - else: - return PhonopyAtoms( - cell=dataset.std_lattice, - scaled_positions=std_positions[perm], - symbols=[atom_data[n][1] for n in std_types[perm]], - magnetic_moments=dataset.std_tensors[perm], - ) - - def get_num_digits(sequence: Sequence, min_length: int = 3) -> int: """Return number of digits of sequence.""" nd = len(str(len(sequence))) diff --git a/src/phelel/velph/utils/plot_eigenvalues.py b/src/phelel/velph/utils/plot_eigenvalues.py index 94d0727..df5764b 100644 --- a/src/phelel/velph/utils/plot_eigenvalues.py +++ b/src/phelel/velph/utils/plot_eigenvalues.py @@ -14,7 +14,7 @@ from phonopy.structure.cells import get_reduced_bases from scipy.spatial import Voronoi -from phelel.velph.cli.utils import get_symmetry_dataset +from phelel.velph.utils.structure import get_symmetry_dataset from phelel.velph.utils.vasp import read_crystal_structure_from_h5 if TYPE_CHECKING: diff --git a/src/phelel/velph/utils/structure.py b/src/phelel/velph/utils/structure.py new file mode 100644 index 0000000..50ad20b --- /dev/null +++ b/src/phelel/velph/utils/structure.py @@ -0,0 +1,137 @@ +"""Utility functions for crystal structure.""" + +from __future__ import annotations + +import numpy as np +import spglib +from numpy.typing import NDArray +from phonopy.interface.vasp import sort_positions_by_symbols +from phonopy.structure.atoms import PhonopyAtoms, get_atomic_data +from phonopy.structure.cells import ( + get_primitive, + get_primitive_matrix_by_centring, + get_reduced_bases, +) +from spglib import SpglibDataset, SpglibMagneticDataset + + +def generate_standardized_cells( + sym_dataset: SpglibDataset | SpglibMagneticDataset, + tolerance: float = 1e-5, +) -> tuple[PhonopyAtoms, PhonopyAtoms, NDArray]: + """Return standardized unit cell and primitive cell.""" + convcell = _get_standardized_unitcell(sym_dataset) + pmat = _get_primitive_matrix_from_dataset(sym_dataset) + if (np.abs(pmat - np.eye(3)) < 1e-8).all(): + primitive = convcell + else: + primitive = get_primitive(convcell, primitive_matrix=pmat, symprec=tolerance) + + return convcell, primitive, pmat + + +def get_primitive_cell( + cell: PhonopyAtoms, + sym_dataset: SpglibDataset | SpglibMagneticDataset, + tolerance: float = 1e-5, +) -> tuple[PhonopyAtoms, NDArray]: + """Return primitive cell and transformation matrix. + + This primitive cell is generated from the input cell without + rigid rotation in contrast to `_get_standardized_unitcell`. + + """ + tmat = sym_dataset.transformation_matrix + pmat = _get_primitive_matrix_from_dataset(sym_dataset) + total_tmat = np.array(np.linalg.inv(tmat) @ pmat, dtype="double", order="C") + + return ( + get_primitive(cell, primitive_matrix=total_tmat, symprec=tolerance), + total_tmat, + ) + + +def get_reduced_cell( + cell: PhonopyAtoms, method="niggli", tolerance: float = 1e-5 +) -> PhonopyAtoms: + """Return a reduced cell of input cell.""" + reduced_lattice = get_reduced_bases(cell.cell, method=method, tolerance=tolerance) + assert reduced_lattice is not None + reduced_positions = cell.scaled_positions @ ( + cell.cell @ np.linalg.inv(reduced_lattice) + ) + reduced_positions = reduced_positions - np.rint(reduced_positions) + reduced_positions[:, :] += np.where(reduced_positions < 0, 1, 0) + reduced_cell = cell.copy() + reduced_cell.cell = reduced_lattice + reduced_cell.scaled_positions = reduced_positions + return reduced_cell + + +def get_symmetry_dataset( + cell: PhonopyAtoms, tolerance: float = 1e-5 +) -> SpglibDataset | SpglibMagneticDataset: + """Return spglib symmetry dataset.""" + if cell.magnetic_moments is None: + dataset = spglib.get_symmetry_dataset(cell.totuple(), symprec=tolerance) + else: + dataset = spglib.get_magnetic_symmetry_dataset( + cell.totuple(), symprec=tolerance + ) + + assert dataset is not None + + return dataset + + +def _get_primitive_matrix_from_dataset( + sym_dataset: SpglibDataset | SpglibMagneticDataset, +) -> NDArray: + """Return primitive matrix from symmetry dataset.""" + if isinstance(sym_dataset, SpglibDataset): + centring = sym_dataset.international[0] + else: + spg_type = spglib.get_spacegroup_type(sym_dataset.hall_number) + assert spg_type is not None + centring = spg_type.international_short[0] + return get_primitive_matrix_by_centring(centring) + + +def _get_standardized_unitcell( + dataset: SpglibDataset | SpglibMagneticDataset, +) -> PhonopyAtoms: + """Return conventional unit cell. + + This conventional unit cell can include rigid rotation with respect to + input unit cell for which symmetry was analized. + + Parameters + ---------- + cell : PhonopyAtoms + Input cell. + dataset : SpgliDataset + Symmetry dataset of spglib. + + Returns + ------- + PhonopyAtoms + Convetional unit cell. + + """ + std_positions = dataset.std_positions + std_types = dataset.std_types + _, _, _, perm = sort_positions_by_symbols(std_types, std_positions) + atom_data = get_atomic_data().atom_data + if isinstance(dataset, SpglibDataset): + return PhonopyAtoms( + cell=dataset.std_lattice, + scaled_positions=std_positions[perm], + symbols=[atom_data[n][1] for n in std_types[perm]], + ) + else: + return PhonopyAtoms( + cell=dataset.std_lattice, + scaled_positions=std_positions[perm], + symbols=[atom_data[n][1] for n in std_types[perm]], + magnetic_moments=dataset.std_tensors[perm], + ) diff --git a/src/phelel/velph/utils/vasp.py b/src/phelel/velph/utils/vasp.py index 442ecb5..0b7473d 100644 --- a/src/phelel/velph/utils/vasp.py +++ b/src/phelel/velph/utils/vasp.py @@ -16,6 +16,105 @@ from spglib import SpgCell +def get_reclat_from_vaspout(fp_vaspout: h5py.File) -> NDArray[np.float64]: + """Return reciprocal space basis vectors. + + Returns + ------- + reclat : np.ndarray + Reciprocal basis vectors in row vectors. + shape=(3, 3) + + """ + # Basis vectors in direct space in column vectors + lattice = np.transpose( + fp_vaspout["input"]["poscar"]["lattice_vectors"][:] # type: ignore + * fp_vaspout["input"]["poscar"]["scale"][()] # type: ignore + ) + # Basis vectors in reciprocal space in row vectors + reclat = 2 * np.pi * np.linalg.inv(lattice) + return reclat + + +def get_bands_data( + kpoint_coords: NDArray, + reclat: NDArray, + nk_per_seg: int, + labels: list[str], +) -> tuple[NDArray, list, list]: + """Measure distances of points from origin along paths and get special points.""" + k_cart = kpoint_coords @ reclat + nk_total = len(kpoint_coords) + n_segments = nk_total // nk_per_seg + assert n_segments * nk_per_seg == nk_total + distances = _get_distances_along_BZ_path(nk_total, n_segments, nk_per_seg, k_cart) + points, labels_at_points = _get_special_points( + labels, distances, n_segments, nk_per_seg, nk_total + ) + + return distances, points, labels_at_points + + +def _get_distances_along_BZ_path( + nk_total: int, n_segments: int, nk_per_seg: int, k_cart: NDArray +) -> NDArray: + """Measure distances of points from origin along paths. + + Returns + ------- + distances : np.ndarray + Distances of points from origin along BZ-paths. + shape=(nk_total,) + + """ + distances = np.zeros(nk_total) + count = 0 + for _ in range(n_segments): + for i_pts in range(nk_per_seg): + # Treatment of jump between equivalent points on BZ boundary + if i_pts == 0: + delta_dist = 0 + else: + delta_dist = np.linalg.norm(k_cart[count, :] - k_cart[count - 1, :]) + distances[count] = distances[count - 1] + delta_dist + count += 1 + return distances + + +def _get_special_points( + labels: list[str], + distances: NDArray, + n_segments: int, + nk_per_seg: int, + nk_total: int, +) -> tuple[list, list]: + """Plot special points at vertical lines and labels.""" + # Left most + points = [] + labels_at_points = [] + + labels_at_points.append(labels[0][0]) + points.append(distances[0]) + + count = 0 + for i_seg in range(n_segments): + for _ in range(nk_per_seg): + count += 1 + if count != nk_total: + points.append(distances[count]) + if labels[i_seg * 2 + 1] == labels[i_seg * 2 + 2]: + labels_at_points.append(labels[i_seg * 2 + 1][0]) + else: + labels_at_points.append( + f"{labels[i_seg * 2 + 1][0]}|{labels[i_seg * 2 + 2][0]}" + ) + + labels_at_points.append(labels[-1][0]) + points.append(distances[-1]) + + return points, labels_at_points + + class VaspKpoints: """KPOINTS file writer. diff --git a/test/velph/cli/el_bands/plot/test_plot.py b/test/velph/cli/el_bands/plot/test_plot.py index 7de1060..37a4965 100644 --- a/test/velph/cli/el_bands/plot/test_plot.py +++ b/test/velph/cli/el_bands/plot/test_plot.py @@ -5,9 +5,10 @@ import h5py import numpy as np import pytest +from numpy.typing import NDArray from phelel.velph.cli.el_bands.plot import _get_bands_data, _get_dos_data -from phelel.velph.cli.utils import get_reclat_from_vaspout +from phelel.velph.utils.vasp import get_reclat_from_vaspout cwd = pathlib.Path(__file__).parent @@ -25,13 +26,14 @@ def test_velph_el_bands_plot_TiNiSn(): 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"], + f_h5py_bands["results/electron_eigenvalues_kpoints_opt"], # type: ignore + f_h5py_bands["input/kpoints_opt"], # type: ignore ) ymin, ymax = 3.575980267703933, 17.575980267703933 - dos, energies, xmax = _get_dos_data( - f_h5py_dos["results/electron_dos_kpoints_opt"], ymin, ymax - ) + f_h5py_dos_results = f_h5py_dos["results/electron_dos_kpoints_opt"] + dos: NDArray = f_h5py_dos_results["dos"][:] # type: ignore + energies: NDArray = f_h5py_dos_results["energies"][:] # type: ignore + xmax = _get_dos_data(dos, energies, ymin, ymax) assert len(distances) == 306 assert pytest.approx(distances[100], 1e-5) == 1.421887803385511 diff --git a/test/velph/cli/init/test_cmd_init.py b/test/velph/cli/init/test_cmd_init.py index a73b26d..092c3d3 100644 --- a/test/velph/cli/init/test_cmd_init.py +++ b/test/velph/cli/init/test_cmd_init.py @@ -39,9 +39,9 @@ VelphFilePaths, VelphInitOptions, VelphInitParams, - get_symmetry_dataset, ) from phelel.velph.templates import default_template_dict +from phelel.velph.utils.structure import get_symmetry_dataset cwd = pathlib.Path(__file__).parent diff --git a/test/velph/cli/test_utils.py b/test/velph/cli/test_utils.py index e7f2ae6..e57fb35 100644 --- a/test/velph/cli/test_utils.py +++ b/test/velph/cli/test_utils.py @@ -9,8 +9,9 @@ from phonopy.interface.phonopy_yaml import read_cell_yaml from phonopy.structure.atoms import PhonopyAtoms -from phelel.velph.cli.utils import get_reduced_cell, get_scheduler_dict +from phelel.velph.cli.utils import get_scheduler_dict from phelel.velph.templates import default_template_dict +from phelel.velph.utils.structure import get_reduced_cell def test_get_scheduler_dict():