Source code for spgrep_modulation.irreps

"""Irreps for phonon."""
from __future__ import annotations

from typing import Literal

import numpy as np
from phonopy.structure.cells import Primitive
from phonopy.structure.symmetry import Symmetry
from spgrep.core import get_spacegroup_irreps_from_primitive_symmetry
from spgrep.representation import project_to_irrep

from spgrep_modulation.utils import NDArrayComplex, NDArrayFloat, NDArrayInt


[docs]def project_eigenmode_representation( eigenmode_representation: NDArrayComplex, primitive: Primitive, primitive_symmetry: Symmetry, primitive_qpoint: NDArrayFloat, method: Literal["Neto", "random"] = "Neto", rtol: float = 1e-5, atol: float = 1e-6, # Too tight tolerance gives wrong basis vectors... ) -> tuple[list[NDArrayComplex], list[NDArrayComplex], NDArrayInt]: """Decompose representation matrices by eigenmodes into irreps. Parameters ---------- eigenmode_representation: array, (order, num_atoms, 3, num_atoms, 3) Representation matrices formed by phonon eigenmodes primitive: phonopy.structure.cells.Primitive phonopy's primitive object primitive_symmetry: phonopy.structure.symmetry.Symmetry phonopy's Symmetry object for primitive cell primitive_qpoint: array, (3, ) q vector in ``primitive``'s dual basis vectors method: str, 'Neto' or 'random' 'Neto': construct irreps from a fixed chain of subgroups of little co-group 'random': construct irreps by numerically diagonalizing a random matrix commute with regular representation rtol: float Relative tolerance for comparing float values atol: float Absolute tolerance to distinguish difference eigenvalues Returns ------- all_basis: list all_basis[i] is list of independent basis vectors with (dim, num_atoms, 3) forming irreps[i]. Note: phase chosen to be consistent with definition of phonopy's dynamical matrix irreps: list of irrep (little_order, dim, dim) mapping_little_group: list of indices for little group. """ rotations = primitive_symmetry.symmetry_operations["rotations"] translations = primitive_symmetry.symmetry_operations["translations"] num_atoms = len(primitive) small_reps, mapping_little_group = get_spacegroup_irreps_from_primitive_symmetry( rotations, translations, primitive_qpoint, method=method, rtol=rtol, atol=atol, ) modified_basis = [] irreps = [] for irrep in small_reps: projected = project_to_irrep( representation=eigenmode_representation[mapping_little_group].reshape( -1, num_atoms * 3, num_atoms * 3 ), irrep=irrep, atol=atol, ) if len(projected) == 0: continue modified_basis.append([basis.reshape(-1, num_atoms, 3) for basis in projected]) irreps.append(irrep) # Apply phase to be consistent with definition of phonopy's dynamical matrix phase = np.exp( -2j * np.pi * np.dot(primitive.scaled_positions, primitive_qpoint) ) # (num_atoms, ) all_basis = [] for list_mb in modified_basis: basis_for_same_irrep = [] for mb in list_mb: eigenvecs = phase[None, :, None] * mb # (dim_irrep, num_atoms, 3) basis_for_same_irrep.append(eigenvecs) all_basis.append(basis_for_same_irrep) return all_basis, irreps, mapping_little_group
[docs]def get_eigenmode_representation( primitive: Primitive, primitive_symmetry: Symmetry, primitive_qpoint: NDArrayFloat, ) -> NDArrayComplex: r"""Compute representation matrix for eigenmodes. .. math:: \Gamma_{\kappa'\mu'; \kappa\mu}^{\mathbf{q}}(g) := \exp \left( -i \mathbf{R}_{g} \mathbf{q} \cdot \mathbf{h}_{g}(\kappa) \right) [\mathbf{R}_{g}]_{\mu'\mu} \delta_{ g\kappa, \kappa' } Parameters ---------- primitive: phonopy.structure.cells.Primitive phonopy's primitive object primitive_symmetry: phonopy.structure.symmetry.Symmetry phonopy's Symmetry object for primitive cell primitive_qpoint: array, (3, ) q vector in ``primitive``'s dual basis vectors Returns ------- rep: array, (order, num_atoms, 3, num_atoms, 3) Representation matrices """ rotations = primitive_symmetry.symmetry_operations["rotations"] translations = primitive_symmetry.symmetry_operations["translations"] order = len(rotations) num_atoms = len(primitive) shifts = np.zeros((order, num_atoms, 3)) for i, (Ri, vi) in enumerate(zip(rotations, translations)): perm_i = primitive_symmetry.atomic_permutations[i] shifts[i] = ( primitive.scaled_positions @ Ri.T + vi[None, :] - primitive.scaled_positions[perm_i] ) perm_rep = np.zeros((order, num_atoms, num_atoms), dtype=np.complex_) for i, Ri in enumerate(rotations): for kappa in range(num_atoms): kappa2 = primitive_symmetry.atomic_permutations[i, kappa] perm_rep[i, kappa2, kappa] = np.exp( -2j * np.pi * np.dot(Ri.T @ primitive_qpoint, shifts[i, kappa]) ) # Rotation matrix in cartesian (order, 3, 3) A = primitive.cell.T # column-wise lattice vectors Ainv = np.linalg.inv(A) rotation_rep = np.array([A @ r @ Ainv for r in rotations], dtype=np.complex_) rep = np.einsum("ipq,iab->ipaqb", perm_rep, rotation_rep, optimize="greedy") return rep