from __future__ import annotations
from argparse import ArgumentParser, RawDescriptionHelpFormatter
import numpy as np
from phonopy.phonon.character_table import character_table as all_character_tables
def build_parser() -> ArgumentParser:
parser = ArgumentParser(
prog="crystod-group --product",
description=(
"Calculate direct products among irreducible representations "
"of a point group."
),
formatter_class=RawDescriptionHelpFormatter,
epilog=(
"Examples:\n"
" crystod-group --product T2g T2g T1u\n --point-group m-3m"
' crystod-group --product E1u E1u --point-group "6/mmm"'
),
)
parser.add_argument(
"--point-group",
"-pg",
required=True,
help="Point group label, e.g. m-3m or 6/mmm.",
)
parser.add_argument(
"--irreps",
"-irreps",
nargs="*",
default=None,
help="Irrep labels to multiply, e.g. T2g T2g T1u.",
)
parser.add_argument(
"--show-irrep-table",
action="store_true",
help="Show the point-group character table.",
)
return parser
def _flatten_irreps(raw_irreps: list[str]) -> list[str]:
irreps: list[str] = []
for item in raw_irreps:
irreps.extend(item.split())
return irreps
def _format_character_value(value) -> str:
scalar = np.asarray(value).item()
if isinstance(scalar, complex):
if abs(scalar.imag) < 1e-10:
scalar = scalar.real
else:
return f"{scalar.real:.4g}{scalar.imag:+.4g}j"
if abs(float(scalar) - round(float(scalar))) < 1e-10:
return str(int(round(float(scalar))))
return f"{float(scalar):.4g}"
def _get_character_table(point_group: str) -> dict:
try:
return all_character_tables[point_group][0]
except KeyError as exc:
available = ", ".join(all_character_tables.keys())
raise SystemExit(
f'ERROR: "{point_group}" is not in the point groups.\n'
f"Choose from: {available}"
) from exc
def format_irrep_table(point_group: str, ct: dict) -> str:
class_names = list(ct["rotation_list"])
class_sizes = [
np.asarray(ct["mapping_table"][class_name]).shape[0]
for class_name in class_names
]
irrep_names = list(ct["character_table"].keys())
header = ["irrep"] + [f"{name}({size})" for name, size in zip(class_names, class_sizes)]
rows = []
for irrep_name in irrep_names:
characters = ct["character_table"][irrep_name]
rows.append(
[irrep_name] + [_format_character_value(value) for value in characters]
)
widths = [len(item) for item in header]
for row in rows:
for idx, item in enumerate(row):
widths[idx] = max(widths[idx], len(item))
lines = []
lines.append(" ".join(item.rjust(widths[idx]) for idx, item in enumerate(header)))
for row in rows:
lines.append(" ".join(item.rjust(widths[idx]) for idx, item in enumerate(row)))
return (
"\n"
"* Point group *\n"
f"{point_group}\n\n"
"* IrRep Table *\n"
"table:\n"
+ "\n".join(lines)
+ "\n"
)
[docs]
def direct_product_character(ct: dict, point_group: str, irreps: list[str]) -> np.ndarray:
"""Characters of the direct product of point-group irreps.
The first step of ``crystod-group --product IRREP... --pg PG``: the
character of a direct product is the product of the characters, class
by class.
Args:
ct: Character table from ``crystod.group.get_character_table``.
point_group: Point-group label (used in the error message only).
irreps: Irrep labels to multiply, e.g. ``["T2g", "T2g", "T1u"]``.
Returns:
The product characters, one per class in the order of
``ct["rotation_list"]``.
Raises:
SystemExit: An irrep label is not in the table; the message lists
the available labels (``ValueError`` when called through
``crystod.group``).
Example:
>>> from crystod import group
>>> ct = group.get_character_table("m-3m")
>>> chi = group.direct_product_character(ct, "m-3m", ["T2g", "T2g"])
>>> chi
array([9., 0., 1., 1., 1., 9., 1., 0., 1., 1.])
>>> {k: n for k, n in group.decompose_representation(ct, chi).items() if n}
{'A1g': 1, 'Eg': 1, 'T1g': 1, 'T2g': 1}
"""
all_irreps = list(ct["character_table"].keys())
irreps_character = []
for irrep in irreps:
if irrep not in all_irreps:
available = ", ".join(all_irreps)
raise SystemExit(
f'ERROR: "{irrep}" is not in irreps of {point_group}.\n'
f"Choose from: {available}"
)
irreps_character.append(np.asarray(ct["character_table"][irrep], dtype=float))
return np.prod(np.stack(irreps_character, axis=0), axis=0)
[docs]
def decompose_representation(ct: dict, reducible_character: np.ndarray) -> dict[str, int]:
"""Reduce a character vector into the irreps of a point group.
The second step of ``crystod-group --product`` with a point group, and
the reduction used by ``--basis``: the reduction formula with the class
sizes taken from the table (same result as ``crystod.group.decompose``,
which takes the class sizes explicitly).
Args:
ct: Character table from ``crystod.group.get_character_table``.
reducible_character: Characters of the reducible representation,
one per class in the order of ``ct["rotation_list"]``.
Returns:
``{irrep label: multiplicity}`` over every irrep of the table (zeros
included).
"""
multiplicities = np.array(
[np.asarray(values).shape[0] for values in ct["mapping_table"].values()],
dtype=float,
)
denominator = float(multiplicities.sum())
results: dict[str, int] = {}
for irrep, character in ct["character_table"].items():
char_array = np.asarray(character, dtype=float)
numerator = np.dot(char_array, multiplicities * reducible_character)
numerator = np.real_if_close(numerator, tol=1000)
if isinstance(numerator, np.ndarray):
numerator = numerator.item()
if isinstance(numerator, complex):
numerator = numerator.real
results[irrep] = round(float(numerator) / denominator)
return results
def format_result(point_group: str, irreps: list[str], results: dict[str, int]) -> str:
formula = "*".join(irreps)
result = " + ".join(f"{value}({key})" for key, value in results.items() if value > 0)
return (
"\n"
"* Point group *\n"
f"{point_group}\n\n"
"* Direct product *\n"
f"{formula}\n\n"
"* Result *\n"
f" {result}\n"
)
def main(argv: list[str] | None = None) -> None:
parser = build_parser()
args = parser.parse_args(argv)
point_group = args.point_group
ct = _get_character_table(point_group)
outputs: list[str] = []
if args.show_irrep_table:
outputs.append(format_irrep_table(point_group, ct).strip("\n"))
irreps = _flatten_irreps(args.irreps or [])
if irreps:
reducible_character = direct_product_character(ct, point_group, irreps)
results = decompose_representation(ct, reducible_character)
outputs.append(format_result(point_group, irreps, results).strip("\n"))
elif not args.show_irrep_table:
parser.error("Specify --irreps and/or --show-irrep-table.")
print("\n\n".join(outputs))
