"""
Star of k: orbit of a k point under the space-group rotations.
Provides a standalone CLI mode and reusable helpers for the modulation mode
(multi-q modulations combine arms of the same star).
"""
from __future__ import annotations
from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser, RawDescriptionHelpFormatter, RawTextHelpFormatter
import numpy as np
from numpy.typing import NDArray
from phonopy.structure.cells import get_primitive_matrix_by_centring
from .operations import get_seitz_symbol
from .runtime_compat import get_little_group
from .spglib_compat import ensure_spglib_compat
ensure_spglib_compat()
from phonopy.interface.calculator import read_crystal_structure
from .vibration_modes import SymmetryOnlyVibrations
class MyHelpFormatter(
RawTextHelpFormatter,
RawDescriptionHelpFormatter,
ArgumentDefaultsHelpFormatter,
):
pass
desc = """
Display the star of k: the set of inequivalent k points generated from a given
k point by the space-group rotations (k' = k R, modulo reciprocal lattice).
# Command Examples:
crystod --star-of-k --poscar 221_PPOSCAR_ScF3 --kpoint 0.5 0.5 0
crystod --star-of-k --poscar 221_PPOSCAR_ScF3 --kpoint M
"""
def build_parser() -> ArgumentParser:
parser = ArgumentParser(description=desc, formatter_class=MyHelpFormatter)
parser.add_argument("--poscar", default="POSCAR", help="POSCAR path.")
parser.add_argument(
"--kpoint",
nargs="+",
required=True,
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.",
)
return parser
def _wrap_to_unit(kpoint: NDArray[np.float64]) -> NDArray[np.float64]:
"""Wrap a k point into [-0.5, 0.5) for display and comparison."""
wrapped = np.remainder(np.asarray(kpoint, dtype=float) + 0.5, 1.0) - 0.5
wrapped[np.abs(wrapped + 0.5) < 1e-8] = 0.5
return wrapped
def _k_equivalent(k_a: NDArray[np.float64], k_b: NDArray[np.float64], atol: float = 1e-8) -> bool:
diff = np.asarray(k_a, dtype=float) - np.asarray(k_b, dtype=float)
return bool((np.abs(diff - np.rint(diff)) < atol).all())
[docs]
def compute_star(
rotations: NDArray[np.int_],
translations: NDArray[np.float64],
kpoint: list[float],
) -> list[dict]:
"""Compute the star of k: the orbit of a k point under the rotations.
Every rotation ``R`` of the space group sends ``k`` to ``k' = k R``;
the distinct images modulo reciprocal-lattice vectors are the arms of
the star, and the rotations sending ``k`` to one arm form a coset of
the little group of ``k``. This is the computation behind
``crystod --star-of-k``; ``crystod-phonon --modulation`` uses it to
combine the arms of a multi-q modulation.
Args:
rotations: Integer rotation matrices of the space group in the
primitive basis, shape ``(n_ops, 3, 3)`` -- for example
``SymmetryAdaptedOrbitalBasis.rotations``.
translations: The matching fractional translations, shape
``(n_ops, 3)``. Accepted for a uniform call signature; the star
depends on the rotations only.
kpoint: Three primitive reciprocal coordinates of ``k``.
Returns:
One dict per arm, in order of first appearance along ``rotations``:
``"kpoint"`` (the arm wrapped into ``[-0.5, 0.5)``),
``"representative_index"`` (index of the first rotation reaching
the arm, the coset representative) and ``"operation_indices"``
(indices of every rotation reaching the arm). With the identity
listed first, as spglib does, the first arm is ``k`` itself and its
``"operation_indices"`` are the little group of ``k``.
Example:
>>> from phonopy.interface.calculator import read_crystal_structure
>>> from crystod import salc
>>> from crystod.examples import example_path
>>> cell, _ = read_crystal_structure(
... str(example_path("221_PPOSCAR_ScF3")), interface_mode="vasp")
>>> basis = salc.SymmetryAdaptedOrbitalBasis(cell=cell)
>>> arms = salc.compute_star(basis.rotations, basis.translations,
... [0.5, 0.5, 0])
>>> [arm["kpoint"].tolist() for arm in arms]
[[0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]]
>>> len(arms[0]["operation_indices"])
16
"""
kpoint = np.asarray(kpoint, dtype=float)
arms: list[dict] = []
for index, rotation in enumerate(rotations):
k_image = kpoint @ rotation
for arm in arms:
if _k_equivalent(k_image, arm["kpoint_raw"]):
arm["operation_indices"].append(index)
break
else:
arms.append(
{
"kpoint_raw": k_image,
"kpoint": _wrap_to_unit(k_image),
"representative_index": index,
"operation_indices": [index],
}
)
for arm in arms:
del arm["kpoint_raw"]
return arms
def print_star_of_k(
rotations: NDArray[np.int_],
translations: NDArray[np.float64],
kpoint: list[float],
seitz_symbols: list[str] | None = None,
indent: str = " ",
) -> list[dict]:
"""Print the star of k and return the computed arms."""
arms = compute_star(rotations, translations, kpoint)
_, _, mapping_little_group = get_little_group(
rotations=rotations,
translations=translations,
kpoint=kpoint,
)
order = len(rotations)
little_order = len(mapping_little_group)
print(f"{indent}|G| = {order}, |G_k| = {little_order}, |star of k| = {len(arms)}")
for line in format_star_lines(arms, seitz_symbols, indent=indent):
print(line)
return arms
def read_poscar_or_exit(poscar_path: str):
"""Read a POSCAR file, exiting with a clear message when it cannot be read."""
import os
if not os.path.isfile(poscar_path):
raise SystemExit(f"ERROR: No POSCAR named {poscar_path}!")
cell, _ = read_crystal_structure(poscar_path, interface_mode="vasp")
if cell is None:
raise SystemExit(f"ERROR: failed to read POSCAR file: '{poscar_path}'")
return cell
def main(argv: list[str] | None = None) -> None:
args = build_parser().parse_args(argv)
cell = read_poscar_or_exit(args.poscar)
structure = SymmetryOnlyVibrations(cell=cell, symprec=args.tolerance)
kpoint_label, kpoint = resolve_kpoint_input(structure, args.kpoint)
dataset = structure.spglib_dataset
print(f"\n * Space group *\n {dataset['international']} ({dataset['number']})\n")
print(f" * k point (primitive) *\n {kpoint_label} {np.round(kpoint, 6).tolist()}\n")
transformation_matrix = get_primitive_matrix_by_centring(dataset["international"][0])
seitz_symbols = [
get_seitz_symbol(rotation, transformation_matrix) for rotation in structure.rotations
]
print(" * Star of k *")
print_star_of_k(
rotations=structure.rotations,
translations=structure.translations,
kpoint=kpoint,
seitz_symbols=seitz_symbols,
)
print("")
if __name__ == "__main__":
main()
