"""
Symmetry-only vibration basis workflow for crystod.
"""
from __future__ import annotations
from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser, RawDescriptionHelpFormatter, RawTextHelpFormatter
from fractions import Fraction
from pathlib import Path
import numpy as np
import spglib
from ase import Atoms
from ase.io import write as ase_write
from numpy.typing import NDArray
from phonopy.structure.cells import get_primitive_matrix_by_centring
from .irreptables_compat import load_irreptables
from .operations import conjugated_little_group_map, find_star_arm, parse_qpoint_token, snap_qpoint
from .runtime_compat import (
SymmetryDatasetAdapter,
get_character,
get_chemical_symbols,
get_little_group,
get_scaled_positions,
)
from .spglib_compat import ensure_spglib_compat
ensure_spglib_compat()
from phonopy.interface.calculator import read_crystal_structure
from phonopy.structure.atoms import PhonopyAtoms
from spgrep.core import get_spacegroup_irreps_from_primitive_symmetry
from spgrep.representation import project_to_irrep
IrrepTable, Irrep = load_irreptables()
class MyHelpFormatter(
RawTextHelpFormatter,
RawDescriptionHelpFormatter,
ArgumentDefaultsHelpFormatter,
):
pass
desc = """
Construct symmetry-allowed vibration basis vectors without phonon force data.
# Command Examples:
crystod-phonon --vibration -c example/test_POSCARs/221_PPOSCAR_ScF3 --qpoint 0.5 0.5 0.5
crystod-phonon --vibration -c example/test_POSCARs/221_PPOSCAR_ScF3 --qpoint R --mode-index 3 --component-index 1 --output POSCAR_vibration
"""
def build_parser() -> ArgumentParser:
parser = ArgumentParser(description=desc, formatter_class=MyHelpFormatter)
parser.add_argument(
"--poscar",
default="POSCAR",
help="POSCAR path.",
)
parser.add_argument(
"--qpoint",
nargs="+",
default=None,
help="Either a high-symmetry label such as GM/X/M/R or three primitive reciprocal coordinates.",
)
parser.add_argument(
"--tolerance",
type=float,
default=1e-5,
help="Symmetry tolerance.",
)
parser.add_argument(
"--list-qpoints",
action="store_true",
help="Only list available high-symmetry q-points and exit.",
)
parser.add_argument(
"--mode-index",
type=int,
default=None,
help="Irrep-grouped mode-space number to inspect (1-based).",
)
parser.add_argument(
"--component-index",
type=int,
default=1,
help="Component number inside the selected degenerate mode space (1-based).",
)
parser.add_argument(
"--amplitude",
type=float,
default=0.3,
help="Amplitude used when writing a displaced structure.",
)
parser.add_argument(
"--output",
default=None,
help="Optional output POSCAR path for the selected mode/component.",
)
parser.add_argument(
"--export-npz",
default=None,
help="Optional .npz path to save positions, displacements, symbols, and lattice.",
)
return parser
class _CoreRepresentation:
def __init__(self, cell: PhonopyAtoms, symprec: float = 1e-5, standardize: bool = True):
if standardize:
primitive_lattice, primitive_pos, primitive_numbers = spglib.standardize_cell(
cell.totuple(),
to_primitive=True,
symprec=symprec,
)
self.primitive_cell = PhonopyAtoms(
numbers=primitive_numbers,
scaled_positions=primitive_pos,
cell=primitive_lattice,
)
print("\n ### Inputed cell was converted into primitive cell. ###")
else:
# Keep the input cell as-is (it must already be primitive). This
# preserves the caller's atom positions so that phase conventions
# stay consistent with an externally built dynamical matrix.
self.primitive_cell = cell
dataset = SymmetryDatasetAdapter(
spglib.get_symmetry_dataset(self.primitive_cell.totuple(), symprec=symprec)
)
self.spglib_dataset = dataset
self.rotations = dataset.rotations
self.translations = dataset.translations
self.symprec = symprec
def _get_isoir_label_list(
self,
qpoint: list[float],
irreps,
mapping_little_group,
) -> list[str] | None:
"""ISO-IR (Miller-Love) labels for spgrep irreps at a non-tabulated q,
or None when the ISO-IR data are unavailable or matching fails.
"""
from .isoir import get_isoir_label_map
self.labels_from_isoir = False
matched = get_isoir_label_map(
self.spglib_dataset["number"],
self.primitive_cell.totuple(),
self.symprec,
qpoint,
[self.rotations[index] for index in mapping_little_group],
[self.translations[index] for index in mapping_little_group],
[get_character(irrep) for irrep in irreps],
)
if matched is None:
return None
label_map, qpoint_name = matched
self.labels_from_isoir = True
self.isoir_qpoint_name = qpoint_name
return [
f"{label_map[index]}({irrep.shape[1]})"
for index, irrep in enumerate(irreps)
]
def get_modified_permutation_rep(
self,
rotation: NDArray[np.int_],
translation: NDArray[np.float64],
kpoint: list[float],
) -> NDArray[np.complex128]:
"""Atom-permutation matrix of one operation at q, with Bloch phases."""
positions = get_scaled_positions(self.primitive_cell)
num_atom = len(positions)
matrix = np.zeros((num_atom, num_atom), dtype=complex)
for i, pos_in in enumerate(positions):
pos_rot = np.dot(rotation, pos_in) + translation
for j, pos_out in enumerate(positions):
diff = pos_rot - pos_out
if (abs(diff - np.rint(diff)) < 1e-5).all():
phase_factor = np.dot(
kpoint,
np.dot(np.linalg.inv(rotation), pos_out - translation) - pos_out,
)
matrix[j, i] = np.exp(2j * np.pi * phase_factor)
return matrix
def get_permutation_reps_at_k(
self,
little_rotations: NDArray[np.int_],
little_translations: NDArray[np.float64],
kpoint: list[float],
) -> NDArray[np.complex128]:
"""Permutation matrices of the little-group operations at q."""
return np.array(
[
self.get_modified_permutation_rep(rotation, translation, kpoint)
for rotation, translation in zip(little_rotations, little_translations)
],
dtype=np.complex128,
)
def get_little_group(self, kpoint: list[float]):
"""Little-group operations of q (see ``runtime_compat.get_little_group``)."""
return get_little_group(
rotations=self.rotations,
translations=self.translations,
kpoint=kpoint,
)
[docs]
class SymmetryOnlyVibrations(_CoreRepresentation):
"""Symmetry-allowed vibration bases of a crystal, without force data.
The engine of ``crystod-phonon --vibration``: at a q point, the
displacement representation of the little group of q (the permutation
representation of the atoms times the Cartesian rotation, with Bloch
phases) is projected onto the spgrep irreps, giving one basis of
symmetry-adapted displacement patterns per irrep occurrence. The spaces
are labeled with ISO-IR irrep names, any component can be written out as
a displaced structure on the commensurate supercell, and
:func:`crystod.phonon.build_symmetry_adapted_modes` uses the same basis to
block-diagonalize a dynamical matrix.
Args:
cell: The crystal structure as a ``phonopy.structure.atoms.PhonopyAtoms``
object, e.g. from ``phonopy.interface.calculator.read_crystal_structure``.
symprec: Symmetry tolerance of the spglib analysis.
standardize: Reduce ``cell`` to the spglib primitive cell first (the
default; a note is printed). ``False`` keeps the input cell as
it is, which must then already be primitive; this preserves the
caller's atom positions so that phase conventions stay consistent
with an externally built dynamical matrix.
Attributes:
primitive_cell: The primitive cell the analysis runs on.
spglib_dataset: Its spglib symmetry dataset (``["number"]``,
``["international"]``, ...).
rotations: Rotation parts of the space-group operations in the
primitive basis, shape ``(n_ops, 3, 3)``.
translations: The corresponding translation parts, shape
``(n_ops, 3)``.
rotations_cartesian: The rotations in Cartesian coordinates.
symprec: The symmetry tolerance.
labels_from_isoir: ``True`` when the last :meth:`get_irrep_labels`
call took its labels from the general ISO-IR k-vector lookup
rather than from the special-point table.
Example:
List the vibration spaces of cubic ScF3 at the R point::
from phonopy.interface.calculator import read_crystal_structure
from crystod import phonon
from crystod.examples import example_path
cell, _ = read_crystal_structure(
example_path("221_PPOSCAR_ScF3"), interface_mode="vasp")
vibrations = phonon.SymmetryOnlyVibrations(cell)
label, qpoint = vibrations.resolve_qpoint(["R"])
irreps, spaces, labels = vibrations.describe_mode_spaces(qpoint)
labels # ['R1+(1)', 'R3+(2)', 'R4+(3)', ...]
[space.shape for space in spaces] # [(1, 12), (2, 12), (3, 12), ...]
"""
def __init__(self, cell: PhonopyAtoms, symprec: float = 1e-5, standardize: bool = True):
super().__init__(cell=cell, symprec=symprec, standardize=standardize)
lattice_t = np.transpose(self.primitive_cell.cell)
lattice_t_inv = np.linalg.inv(lattice_t)
self.rotations_cartesian = np.array(
[lattice_t @ rotation @ lattice_t_inv for rotation in self.rotations],
dtype=np.complex128,
)
[docs]
def get_high_symmetry_qpoints(self) -> dict[str, list[float]]:
"""High-symmetry q points of the primitive cell, from seekpath.
``crystod-phonon --vibration --list-qpoints`` prints this map. The
coordinates are expressed in the reciprocal basis of *this* object's
(spglib) primitive cell: seekpath's own primitive cell can differ
from it by an integer change of basis (base-centred monoclinic cells,
for instance), and the tabulated coordinates are transformed
accordingly. A warning is issued only when the two cells are not
related by such a change of basis, in which case the coordinates are
returned as seekpath gives them.
Returns:
Dict mapping seekpath labels (``"GAMMA"``, ``"R"``, ``"X"``, ...)
to fractional coordinates in the primitive reciprocal basis.
"""
import seekpath
import warnings
structure = (
self.primitive_cell.cell,
self.primitive_cell.scaled_positions,
self.primitive_cell.numbers,
)
path_data = seekpath.get_path(structure, symprec=1e-5)
point_coords = path_data["point_coords"]
own_lattice = np.asarray(self.primitive_cell.cell, dtype=float)
seekpath_lattice = np.asarray(path_data["primitive_lattice"], dtype=float)
if not np.allclose(own_lattice, seekpath_lattice, atol=1e-4):
# rows of the seekpath cell as integer combinations of our rows:
# P_seek = M P_own; a k point with fractional coordinates k_seek
# in seekpath's reciprocal basis is k_own = k_seek inv(M)^T here
change = seekpath_lattice @ np.linalg.inv(own_lattice)
rounded = np.rint(change)
if (np.allclose(change, rounded, atol=1e-6)
and abs(round(np.linalg.det(rounded))) == 1):
to_own = np.linalg.inv(rounded).T
point_coords = {
label: (np.asarray(coords, dtype=float) @ to_own).tolist()
for label, coords in point_coords.items()
}
else:
warnings.warn(
"The primitive cell from seekpath does not match the spglib "
"primitive cell. The q-point coordinates might need a basis "
"transformation.",
stacklevel=2,
)
return point_coords
[docs]
def resolve_qpoint(self, raw_qpoint: list[str]) -> tuple[str, list[float]]:
"""Resolve ``--qpoint`` tokens into a label and coordinates.
One token is a seekpath label (``GM``, ``G`` and the Greek capital
gamma are accepted for ``GAMMA``); three tokens are coordinates in
the primitive reciprocal basis, fractions such as ``1/3`` allowed.
Coordinates are labeled with the special point they coincide with,
else with the name of the star arm the space-group rotations map them
onto, else with the ISO-IR k-vector type of q, else ``"custom"``.
Args:
raw_qpoint: The tokens, one label or three coordinate strings.
Returns:
``(label, qpoint)`` with ``qpoint`` a list of three floats.
Raises:
ValueError: For an unknown label, or a token count other than one
or three.
"""
qpoint_map = self.get_high_symmetry_qpoints()
alias_map = {
"GM": "GAMMA",
"G": "GAMMA",
"Γ": "GAMMA",
}
if len(raw_qpoint) == 1:
requested = raw_qpoint[0].strip().upper()
# the aliases stand in for GAMMA only where the cell has no
# special point of that name (body-centred tetragonal cells have
# a genuine G)
if requested not in qpoint_map:
requested = alias_map.get(requested, requested)
if requested in qpoint_map:
return requested, list(qpoint_map[requested])
available = ", ".join(sorted(qpoint_map))
raise ValueError(
f"Unknown q-point label '{raw_qpoint[0]}'. Available labels: {available}"
)
if len(raw_qpoint) != 3:
raise ValueError("--qpoint must be either one label or three coordinates.")
qpoint = [parse_qpoint_token(value) for value in raw_qpoint]
matched_label = None
for label, coords in qpoint_map.items():
if np.allclose(qpoint, coords, atol=1e-8):
matched_label = label
break
if matched_label is None:
# q may be a non-tabulated arm of a special-point star: label it
# with the name of the arm the space-group rotations map it onto.
arm = find_star_arm(qpoint, self.rotations, list(qpoint_map.values()))
if arm is not None:
for label, coords in qpoint_map.items():
if np.allclose(arm[1], coords, atol=1e-8):
matched_label = label
break
if matched_label is None:
# non-special q: fall back to the ISO-IR k-vector type label
from .isoir import get_isoir_kpoint_name
matched_label = get_isoir_kpoint_name(
self.spglib_dataset["number"], self.primitive_cell.totuple(),
self.symprec, qpoint,
)
return matched_label or "custom", qpoint
[docs]
def get_vibration_rep(self, kpoint: list[float]):
"""Displacement representation of the little group of q.
Args:
kpoint: Fractional coordinates of q in the primitive reciprocal
basis.
Returns:
``(irreps, vibration_rep, mapping_little_group)``: the spgrep
irreps of the little group of q; the representation matrices of
the little-group operations on the ``3 * n_atoms`` displacement
space, shape ``(n_little, 3 * n_atoms, 3 * n_atoms)``; and the
indices of the little-group operations within ``rotations``.
"""
irreps, mapping_little_group = get_spacegroup_irreps_from_primitive_symmetry(
rotations=self.rotations,
translations=self.translations,
kpoint=kpoint,
)
little_rotations = self.rotations[mapping_little_group]
little_translations = self.translations[mapping_little_group]
permutation_matrices = self.get_permutation_reps_at_k(
little_rotations=little_rotations,
little_translations=little_translations,
kpoint=kpoint,
)
cartesian_rep = self.rotations_cartesian[mapping_little_group]
vibration_rep = np.array(
[
np.kron(permutation_matrix, cartesian_rotation)
for permutation_matrix, cartesian_rotation in zip(permutation_matrices, cartesian_rep)
],
dtype=np.complex128,
)
return irreps, vibration_rep, mapping_little_group
[docs]
def get_vibration_basis(
self,
irreps,
vibration_rep,
irrep_labels: list[str] | None = None,
) -> tuple[list[NDArray[np.complex128]], list[str]]:
"""Project the displacement representation onto each irrep.
Args:
irreps: The spgrep irreps from :meth:`get_vibration_rep`.
vibration_rep: The representation matrices from
:meth:`get_vibration_rep`.
irrep_labels: One label per irrep, e.g. from
:meth:`get_irrep_labels`; generic ``irrep_N(dim)`` labels are
used when omitted.
Returns:
``(basis_vectors, basis_labels)``: one ``(dim, 3 * n_atoms)`` array
per occurrence of an irrep in the displacement representation
(rows are the symmetry-adapted displacement patterns, as bras),
and the label of each space.
"""
basis_vectors: list[NDArray[np.complex128]] = []
basis_labels: list[str] = []
fallback_labels = irrep_labels or [f"irrep_{index + 1}({irrep.shape[1]})" for index, irrep in enumerate(irreps)]
for irrep, irrep_label in zip(irreps, fallback_labels):
projected_spaces = project_to_irrep(vibration_rep, irrep)
basis_vectors.extend(projected_spaces)
basis_labels.extend([irrep_label] * len(projected_spaces))
return basis_vectors, basis_labels
def _get_irt_irreps_at_q(self, qpoint: list[float], irt_table, prim_mat) -> list[Irrep]:
irreps_at_q = []
prim_inv = np.linalg.inv(prim_mat)
conventional_q = np.array(qpoint) @ prim_inv
for irrep_at_q in irt_table.irreps:
if np.allclose(irrep_at_q.k, conventional_q):
irreps_at_q.append(irrep_at_q)
return irreps_at_q
def _get_mapping_to_irt(
self,
irt_little_rotations: NDArray[np.int_],
found_little_rotations: NDArray[np.int_],
prim_mat: NDArray[np.float64],
) -> list[int]:
conventional_little_rotations = prim_mat @ found_little_rotations @ np.linalg.inv(prim_mat)
mapping_to_irt = []
for irt_rotation in irt_little_rotations:
for index, rotation in enumerate(conventional_little_rotations):
if np.allclose(irt_rotation, rotation):
mapping_to_irt.append(index)
break
return mapping_to_irt
[docs]
def get_irrep_labels(
self,
qpoint: list[float],
irreps,
mapping_little_group: NDArray[np.int_],
) -> list[str]:
"""ISO-IR labels of the spgrep irreps at q.
The characters of each spgrep irrep are matched against the ISO-IR
table of the space group: directly at a tabulated special point, by
conjugation onto the tabulated arm for another arm of its star, and
through the general ISO-IR k-vector lookup (Miller-Love labels) for a
symmetry line, plane or generic q. An irrep no table matches keeps its
generic ``irrep_N(dim)`` label.
Args:
qpoint: Fractional coordinates of q in the primitive reciprocal
basis.
irreps: The spgrep irreps from :meth:`get_vibration_rep`.
mapping_little_group: The little-group indices from
:meth:`get_vibration_rep`.
Returns:
One label per irrep, e.g. ``"R4+(3)"`` (the irrep name with its
dimension).
"""
generic_labels = [f"irrep_{index + 1}({irrep.shape[1]})" for index, irrep in enumerate(irreps)]
self.labels_from_isoir = False
try:
irt_table = IrrepTable(self.spglib_dataset["number"], spinor=False)
except Exception:
return (
self._get_isoir_label_list(qpoint, irreps, mapping_little_group)
or generic_labels
)
prim_mat = get_primitive_matrix_by_centring(self.spglib_dataset["international"][0])
irt_irreps = self._get_irt_irreps_at_q(qpoint, irt_table, prim_mat)
conjugated = None
if not irt_irreps:
# q may be a non-tabulated arm of a special-point star: map it onto
# the tabulated arm and transport the characters by conjugation.
special_points: list[list[float]] = []
for irrep_at_q in irt_table.irreps:
primitive_q = snap_qpoint(np.array(irrep_at_q.k) @ prim_mat)
if primitive_q not in special_points:
special_points.append(primitive_q)
arm = find_star_arm(qpoint, self.rotations, special_points)
if arm is not None:
candidate_irreps = self._get_irt_irreps_at_q(arm[1], irt_table, prim_mat)
transported = conjugated_little_group_map(
self.rotations, self.translations, arm[0], arm[1], mapping_little_group
)
if candidate_irreps and transported is not None:
irt_irreps = candidate_irreps
conjugated = transported
if not irt_irreps:
# Not tabulated as a special point (e.g. a symmetry line/plane or
# generic q): fall back to the general ISO-IR (ISOTROPY) k-vector
# lookup, which covers every k-vector type. Labels then follow
# the Miller-Love convention.
return (
self._get_isoir_label_list(qpoint, irreps, mapping_little_group)
or generic_labels
)
irt_little_rotations = np.array(
[irt_table.symmetries[index - 1].R for index in irt_irreps[0].characters.keys()]
)
found_little_rotations = self.rotations[mapping_little_group]
if conjugated is None:
mapping_to_irt = self._get_mapping_to_irt(irt_little_rotations, found_little_rotations, prim_mat)
character_phases = np.ones(len(irt_little_rotations), dtype=complex)
else:
# mapping_to_irt[m] = position (within the little group of q) of the
# operation h whose conjugate g^-1 h g is the m-th tabulated
# operation; the transported character picks up the Bloch phase.
conj_indices, conj_phases = conjugated
prim_mat_inv = np.linalg.inv(prim_mat)
mapping_to_irt = []
phases: list[complex] = []
for irt_rotation in irt_little_rotations:
rotation_prim = np.rint(prim_mat_inv @ irt_rotation @ prim_mat).astype(int)
table_op_index = None
for j, rotation in enumerate(self.rotations):
if (rotation == rotation_prim).all():
table_op_index = j
break
if table_op_index is None or table_op_index not in conj_indices:
break
position = conj_indices.index(table_op_index)
mapping_to_irt.append(position)
phases.append(conj_phases[position])
character_phases = np.array(phases, dtype=complex)
if len(mapping_to_irt) != len(irt_little_rotations):
return (
self._get_isoir_label_list(qpoint, irreps, mapping_little_group)
or generic_labels
)
resolved_labels: list[str] = []
used_irt_labels: set[str] = set()
for generic_label, irrep in zip(generic_labels, irreps):
spgrep_character = np.array(get_character(irrep), dtype=complex)[mapping_to_irt]
best_label = generic_label
best_overlap = -1.0
for irt_irrep in irt_irreps:
irt_label = f"{irt_irrep.name}({irt_irrep.dim})"
irt_character = (
np.array(list(irt_irrep.characters.values()), dtype=complex) * character_phases
)
overlap = np.abs(
np.dot(spgrep_character, np.conjugate(irt_character)) / irt_irrep.nsym
)
if overlap > best_overlap:
best_label = irt_label
best_overlap = float(overlap)
if best_overlap < 0.9:
best_label = generic_label
elif best_label in used_irt_labels:
best_label = f"{best_label} [{generic_label}]"
used_irt_labels.add(best_label)
resolved_labels.append(best_label)
return resolved_labels
[docs]
def describe_mode_spaces(
self,
qpoint: list[float],
) -> tuple[object, list[NDArray[np.complex128]], list[str]]:
"""Irreps, projected vibration spaces and their labels at q.
The one-call form of :meth:`get_vibration_rep`,
:meth:`get_irrep_labels` and :meth:`get_vibration_basis`, as
``crystod-phonon --vibration`` prints them (one "Mode Space" line per
irrep occurrence, with its label and dimension).
Args:
qpoint: Fractional coordinates of q in the primitive reciprocal
basis.
Returns:
``(irreps, basis_spaces, basis_space_labels)``: the spgrep irreps,
one ``(dim, 3 * n_atoms)`` array per irrep occurrence, and the
ISO-IR label of each space. The mode-space numbers of the command
are 1-based positions in ``basis_spaces``.
"""
irreps, vibration_rep, mapping_little_group = self.get_vibration_rep(qpoint)
irrep_labels = self.get_irrep_labels(qpoint, irreps, mapping_little_group)
basis_spaces, basis_space_labels = self.get_vibration_basis(irreps, vibration_rep, irrep_labels)
return irreps, basis_spaces, basis_space_labels
[docs]
def get_supercell_size(self, qpoint: list[float]) -> tuple[int, int, int]:
"""Supercell multiplicities along a, b, c commensurate with q.
Args:
qpoint: Fractional coordinates of q in the primitive reciprocal
basis.
Returns:
``(n1, n2, n3)``: 1 for a zero component, else the denominator of
the component (limited to 6).
"""
sizes = []
for component in qpoint:
if abs(component) < 1e-10:
sizes.append(1)
else:
sizes.append(Fraction(float(component)).limit_denominator(6).denominator)
return tuple(sizes)
[docs]
def get_supercell_displacements(
self,
qpoint: list[float],
mode_vector: NDArray[np.complex128],
supercell_size: tuple[int, int, int],
):
"""Displacement pattern of one basis vector on a supercell.
Atom j of the primitive cell at lattice translation R is displaced by
``Re(mode_j * exp(2 pi i q . R))`` (unit amplitude).
Args:
qpoint: Fractional coordinates of q in the primitive reciprocal
basis.
mode_vector: One row of a projected space, ``3 * n_atoms``
complex components.
supercell_size: ``(n1, n2, n3)`` multiplicities, e.g. from
:meth:`get_supercell_size`.
Returns:
``(positions, displacements, symbols, supercell_lattice)``:
Cartesian positions and displacements of the supercell atoms,
shape ``(n1 * n2 * n3 * n_atoms, 3)``, their chemical symbols, and
the supercell lattice vectors as rows.
"""
primitive = self.primitive_cell
n_atoms = len(primitive.scaled_positions)
lattice = primitive.cell
frac_pos = primitive.scaled_positions
symbols_prim = get_chemical_symbols(primitive)
mode = mode_vector.reshape((-1, 3))
n1, n2, n3 = supercell_size
all_positions = []
all_displacements = []
all_symbols = []
for i1 in range(n1):
for i2 in range(n2):
for i3 in range(n3):
translation_frac = np.array([i1, i2, i3])
phase = np.exp(2j * np.pi * np.dot(qpoint, translation_frac))
for atom_index in range(n_atoms):
pos_frac = frac_pos[atom_index] + translation_frac
pos_cart = pos_frac @ lattice
all_positions.append(pos_cart)
displacement = np.real(mode[atom_index] * phase)
all_displacements.append(displacement)
all_symbols.append(symbols_prim[atom_index])
supercell_lattice = lattice.copy()
supercell_lattice[0] *= n1
supercell_lattice[1] *= n2
supercell_lattice[2] *= n3
return (
np.array(all_positions),
np.array(all_displacements),
all_symbols,
supercell_lattice,
)
[docs]
def write_displaced_structure(
self,
positions: NDArray[np.float64],
displacements: NDArray[np.float64],
symbols: list[str],
supercell_lattice: NDArray[np.float64],
amplitude: float,
output_path: str,
) -> None:
"""Write ``positions + amplitude * displacements`` as a POSCAR.
Args:
positions: Cartesian positions from
:meth:`get_supercell_displacements`.
displacements: The unit-amplitude displacements from the same
call.
symbols: Chemical symbols of the atoms.
supercell_lattice: Supercell lattice vectors as rows.
amplitude: Displacement amplitude in Angstroms.
output_path: Output file path (VASP format, direct coordinates).
"""
atoms = Atoms(
symbols=symbols,
positions=positions + amplitude * displacements,
cell=supercell_lattice,
pbc=True,
)
ase_write(output_path, atoms, format="vasp", direct=True)
def _print_high_symmetry_qpoints(qpoints: dict[str, list[float]]) -> None:
print("Available high-symmetry q-points:")
for label, coords in qpoints.items():
print(f" {label:8s} {coords}")
def _print_mode_spaces(basis_spaces: list[NDArray[np.complex128]], irrep_labels: list[str]) -> None:
print("Irrep-grouped vibration spaces:")
for mode_index, (space, irrep_label) in enumerate(zip(basis_spaces, irrep_labels), start=1):
dim = space.shape[0]
print(
f" Mode Space {mode_index:2d}: irrep = {irrep_label}, dimension = {dim}, "
f"component numbers = 1..{dim}"
)
def main(argv: list[str] | None = None) -> None:
from .star_of_k import read_poscar_or_exit
args = build_parser().parse_args(argv)
cell = read_poscar_or_exit(args.poscar)
vibrations = SymmetryOnlyVibrations(cell=cell, symprec=args.tolerance)
qpoint_map = vibrations.get_high_symmetry_qpoints()
_print_high_symmetry_qpoints(qpoint_map)
if args.list_qpoints:
return
if not args.qpoint:
raise ValueError("--qpoint is required unless --list-qpoints is used.")
qpoint_label, qpoint = vibrations.resolve_qpoint(args.qpoint)
print(f"\nSelected q-point: {qpoint_label} = {qpoint}")
irreps, basis_spaces, irrep_labels = vibrations.describe_mode_spaces(qpoint)
print(f"Number of irrep-grouped vibration spaces: {len(basis_spaces)}")
_print_mode_spaces(basis_spaces, irrep_labels)
if args.mode_index is None:
print(
"\nUse --mode-index and optionally --component-index to inspect a specific basis vector."
)
return
if args.mode_index < 1 or args.mode_index > len(basis_spaces):
raise SystemExit(
f"ERROR: mode-space number {args.mode_index} is out of range "
f"[1, {len(basis_spaces)}] (numbering is 1-based)."
)
selected_space = basis_spaces[args.mode_index - 1]
if args.component_index < 1 or args.component_index > selected_space.shape[0]:
raise SystemExit(
f"ERROR: component number {args.component_index} is out of range "
f"[1, {selected_space.shape[0]}] (numbering is 1-based)."
)
mode_vector = selected_space[args.component_index - 1]
supercell_size = vibrations.get_supercell_size(qpoint)
print(f"\nSelected mode space: {args.mode_index}")
print(f"Selected irrep : {irrep_labels[args.mode_index - 1]}")
print(f"Selected component : {args.component_index}")
print(f"Commensurate supercell size: {supercell_size}")
positions, displacements, symbols, supercell_lattice = vibrations.get_supercell_displacements(
qpoint=qpoint,
mode_vector=mode_vector,
supercell_size=supercell_size,
)
norms = np.linalg.norm(displacements, axis=1)
print(f"Supercell atom count: {len(symbols)}")
print(f"Displacement norm range (unit amplitude): min={norms.min():.6f}, max={norms.max():.6f}")
print("First 5 displacement vectors:")
for index in range(min(5, len(displacements))):
print(
f" {index:2d} {symbols[index]:2s} "
f"pos={np.round(positions[index], 6).tolist()} "
f"disp={np.round(displacements[index], 6).tolist()}"
)
if args.export_npz:
np.savez(
args.export_npz,
positions=positions,
displacements=displacements,
symbols=np.array(symbols, dtype=object),
supercell_lattice=supercell_lattice,
qpoint=np.array(qpoint, dtype=float),
qpoint_label=np.array(qpoint_label, dtype=object),
mode_index=np.array(args.mode_index),
component_index=np.array(args.component_index),
irrep_labels=np.array(irrep_labels, dtype=object),
selected_irrep_label=np.array(irrep_labels[args.mode_index - 1], dtype=object),
mode_space_dimensions=np.array([space.shape[0] for space in basis_spaces], dtype=int),
selected_mode_dimension=np.array(selected_space.shape[0], dtype=int),
amplitude=np.array(args.amplitude, dtype=float),
)
print(f"Saved mode data to: {args.export_npz}")
if args.output:
vibrations.write_displaced_structure(
positions=positions,
displacements=displacements,
symbols=symbols,
supercell_lattice=supercell_lattice,
amplitude=args.amplitude,
output_path=args.output,
)
print(f"Saved displaced structure to: {args.output}")
if __name__ == "__main__":
main()
