# find_point_group.py (refactored to use tkpointgroup.py)
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
概要:
    XYZ 分子から点群を推定し、対称操作を直交化、整列、ラベル付けして表示します。
詳細説明:
    このスクリプトは、XYZ分子ファイルから点群を推定し、検出された対称操作を分類・ラベル付けするツールです。
    質量中心への平行移動、慣性主軸への整列などの前処理オプションも提供します。
    pymatgen を使用してXYZファイルを読み込み、PointGroupAnalyzer で点群記号と対称操作を特定します。
    対称操作の行列はSVDにより直交化され、tkpointgroup ライブラリの関数を用いて代表値にスナップされます。
    Cnv (n=3,4,6), C2v, Dn (n=3,4,6), T, Td, Th, O, Oh などの点群ファミリーに対応するカスタムラベリング関数により、
    各対称操作に詳細なSchoenflies記号のラベルが付与されます。
    Cnv群ではσvとσdの区別、C2v群ではσv(xz)とσv'(yz)の区別、Dn群ではCn^k(z)とC2'(x系)の分類が行われます。
    必要となるライブラリは pymatgen です。pip install pymatgen でインストールしてください。
    tkpointgroup.py はこのスクリプトと同じディレクトリ、またはPYTHONPATHに存在する必要があります。
関連リンク:
    find_point_group_usage
"""

# コマンドライン実行例はDocstringではなく、Pythonのコメントとして保持します。
# python find_point_group.py --xyz h2o.xyz --center --align --dump-ops --write-aligned h2o_aligned.xyz
# python find_point_group.py --xyz nh3.xyz --center --align --symprec 2e-3 --angle-tol 8
# python find_point_group.py --xyz ch4.xyz --center --align --dump-ops
# python find_point_group.py --xyz h2o.xyz --center --align 
# python find_point_group.py --xyz SF6.xyz --center --align 
# python find_point_group.py --xyz C6H6.xyz --center --align 
# 推奨：オプションはそのまま（eigen-tol=1e-2, matrix-tol=1e-1）

import argparse
import numpy as np
from pathlib import Path
from typing import List, Tuple, Optional

from pymatgen.core import Molecule
from pymatgen.core.periodic_table import Element
from pymatgen.symmetry.analyzer import PointGroupAnalyzer

import tkpointgroup as pg  # ← 追加：共通ライブラリ

# ====== 数値ユーティリティ（pg を用いつつ、このスクリプト特有の処理を保持） ======

def _orth_preserve_det(R: np.ndarray) -> np.ndarray:
    """
    概要:
        行列を直交化し、detの符号を保ちつつ代表的な値にスナップします。
    詳細説明:
        SVDを用いて行列Rを直交行列に変換し、元の行列式detの符号を維持するように調整します。
        最終的には tkpointgroup.snap_matrix を使用して要素を代表値にスナップします。
    引数:
        :param R: 直交化する行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: 直交化され、代表値にスナップされた行列。
        :rtype: numpy.ndarray
    """
    U, _, Vt = np.linalg.svd(R)
    R_ = U @ Vt
    if np.sign(np.linalg.det(R_)) != np.sign(np.linalg.det(R)):
        U[:, -1] *= -1
        R_ = U @ Vt
    return pg.snap_matrix(R_)  # 最終スナップはライブラリを使用

def _vec_to_unit(v: np.ndarray) -> np.ndarray:
    """
    概要:
        ベクトルを単位ベクトルに正規化します。
    詳細説明:
        ベクトルのノルムが非常に小さい場合は、ゼロベクトルを返します。
    引数:
        :param v: 正規化するベクトル。
        :type v: numpy.ndarray
    戻り値:
        :returns: 単位ベクトル。
        :rtype: numpy.ndarray
    """
    n = np.linalg.norm(v)
    return v / n if n > 1e-15 else v

def rot_axis_angle(R: np.ndarray) -> Tuple[Optional[np.ndarray], float]:
    """
    概要:
        回転行列 R の回転軸（Noneのときは不定）と回転角(0..π)を返します。
    詳細説明:
        回転行列から固有値と固有ベクトルを計算し、固有値1に対応する固有ベクトルを回転軸とします。
        回転角はトレースから計算します。回転軸が不定の場合はNoneを返します。
    引数:
        :param R: 回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: 回転軸（Noneの場合あり）と回転角（ラジアン、0からπの範囲）のタプル。
        :rtype: tuple
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-10):
        return None, 0.0
    vals, vecs = np.linalg.eig(R)
    idx = np.where(np.isclose(vals.real, 1.0, atol=1e-6))[0]
    if idx.size == 0:
        th = float(np.arccos(np.clip((np.trace(R)-1)/2, -1, 1)))
        return None, th
    axis = vecs[:, idx[0]].real
    axis /= (np.linalg.norm(axis) + 1e-15)
    th = float(np.arccos(np.clip((np.trace(R)-1)/2, -1, 1)))
    return axis, th

def mirror_normal(R: np.ndarray) -> Optional[np.ndarray]:
    """
    概要:
        鏡映の法線ベクトルを返します（鏡映でない場合は None）。
    詳細説明:
        行列式が負でかつ二乗して単位行列となる行列を鏡映と判定します。
        鏡映の場合、固有値-1に対応する固有ベクトルを法線ベクトルとします。
    引数:
        :param R: 検査する行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: 鏡映の法線ベクトル、またはNone。
        :rtype: Optional[numpy.ndarray]
    """
    R = _orth_preserve_det(R)
    if not (np.linalg.det(R) < 0 and np.allclose(R @ R, np.eye(3), atol=1e-6)):
        return None
    vals, vecs = np.linalg.eig(R)
    idx = np.where(np.isclose(vals.real, -1.0, atol=1e-6))[0]
    if idx.size == 0:
        return None
    n = vecs[:, idx[0]].real
    n /= (np.linalg.norm(n) + 1e-15)
    return n

def pretty_matrix(R: np.ndarray) -> str:
    """
    概要:
        行列を整形して文字列として返します。
    詳細説明:
        微小な要素を0として丸め、指定された書式で数値文字列を生成します。
    引数:
        :param R: 整形する行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: 整形された行列の文字列表現。
        :rtype: str
    """
    Rt = np.where(np.abs(R) < 1e-10, 0.0, R)
    Rt = np.round(Rt, 6)
    return np.array2string(Rt, formatter={'float_kind': lambda x: f"{x:9.6f}"})

# 代表方向（規格化済み）
U100 = [_vec_to_unit(np.array([sx,0,0])) for sx in (1,-1)] + \
       [_vec_to_unit(np.array([0,sy,0])) for sy in (1,-1)] + \
       [_vec_to_unit(np.array([0,0,sz])) for sz in (1,-1)]

U110: List[np.ndarray] = []
for sx in (1,-1):
    for sy in (1,-1):
        U110 += [_vec_to_unit(np.array([sx, sy, 0]))]
for sy in (1,-1):
    for sz in (1,-1):
        U110 += [_vec_to_unit(np.array([0, sy, sz]))]
for sx in (1,-1):
    for sz in (1,-1):
        U110 += [_vec_to_unit(np.array([sx, 0, sz]))]

U111 = [_vec_to_unit(np.array([sx, sy, sz]))
        for sx in (1,-1) for sy in (1,-1) for sz in (1,-1)]

def align_score(axis: np.ndarray, fam: List[np.ndarray]) -> float:
    """
    概要:
        特定の方向群に対する軸の整列スコアを計算します。
    詳細説明:
        与えられた軸が方向群のいずれかのベクトルとどれだけ平行に近いかを示すスコアを計算します。
        最大の内積の絶対値をスコアとします。
    引数:
        :param axis: 評価する軸。
        :type axis: numpy.ndarray
        :param fam: 比較する方向ベクトル群。
        :type fam: List[numpy.ndarray]
    戻り値:
        :returns: 整列スコア (0.0から1.0)。
        :rtype: float
    """
    axis = _vec_to_unit(axis)
    return max(abs(float(np.dot(axis, d))) for d in fam)

# ====== 幾何前処理（共通関数はこのスクリプトで保持） ======

def center_of_mass(mol: Molecule) -> np.ndarray:
    """
    概要:
        分子の質量中心を計算します。
    引数:
        :param mol: pymatgen.core.Molecule オブジェクト。
        :type mol: pymatgen.core.Molecule
    戻り値:
        :returns: 分子の質量中心座標。
        :rtype: numpy.ndarray
    """
    masses = np.array([float(Element(s.species_string).atomic_mass) for s in mol.sites])
    coords = np.array([s.coords for s in mol.sites])
    return (masses[:, None] * coords).sum(axis=0) / masses.sum()

def inertia_align_matrix(mol: Molecule) -> np.ndarray:
    """
    概要:
        分子の慣性主軸に整列するための回転行列を計算します。
    詳細説明:
        質量中心を計算し、それに基づいて慣性テンソルを構築します。
        慣性テンソルの固有ベクトルが慣性主軸となり、これらを基底とする回転行列を生成します。
        行列式が負になる場合は、主軸の一つを反転して右手系を保ちます。
    引数:
        :param mol: pymatgen.core.Molecule オブジェクト。
        :type mol: pymatgen.core.Molecule
    戻り値:
        :returns: 慣性主軸に整列させるための回転行列。
        :rtype: numpy.ndarray
    """
    masses = np.array([float(Element(s.species_string).atomic_mass) for s in mol.sites])
    coords = np.array([s.coords for s in mol.sites])
    com = center_of_mass(mol)
    r = coords - com
    I = np.zeros((3, 3))
    for m, v in zip(masses, r):
        x, y, z = v
        I += m * np.array([[y*y+z*z, -x*y,    -x*z],
                           [-x*y,     x*x+z*z, -y*z ],
                           [-x*z,     -y*z,    x*x+y*y]])
    _, vecs = np.linalg.eigh(I)
    R = vecs.copy()
    if np.linalg.det(R) < 0:
        R[:, 0] *= -1
    return R

def apply_rigid_transform(mol: Molecule, R: np.ndarray, t: np.ndarray) -> Molecule:
    """
    概要:
        分子に剛体変換（回転と平行移動）を適用します。
    詳細説明:
        分子の原子座標に平行移動 t を適用し、その後に回転行列 R を適用します。
    引数:
        :param mol: 変換を適用する pymatgen.core.Molecule オブジェクト。
        :type mol: pymatgen.core.Molecule
        :param R: 回転行列。
        :type R: numpy.ndarray
        :param t: 平行移動ベクトル。
        :type t: numpy.ndarray
    戻り値:
        :returns: 変換が適用された新しい Molecule オブジェクト。
        :rtype: pymatgen.core.Molecule
    """
    coords = np.array([s.coords for s in mol.sites])
    new = (coords - t) @ R.T
    species = [s.species for s in mol.sites]
    return Molecule(species, new)

def write_xyz(mol: Molecule, path: Path, comment: str = ""):
    """
    概要:
        Molecule オブジェクトをXYZファイルとして書き込みます。
    引数:
        :param mol: 書き込む Molecule オブジェクト。
        :type mol: pymatgen.core.Molecule
        :param path: 出力ファイルのパス。
        :type path: pathlib.Path
        :param comment: XYZファイルの2行目に書き込むコメント。
        :type comment: str
    戻り値:
        なし
    """
    with open(path, "w", encoding="utf-8") as f:
        f.write(f"{len(mol)}\n{comment}\n")
        for site in mol.sites:
            el = site.species_string
            x, y, z = site.coords
            f.write(f"{el:2s} {x: .6f} {y: .6f} {z: .6f}\n")

# ====== 点群推定 ======

def guess_point_group(mol: Molecule, tolerance: float = 0.03, eigen_tolerance: float = 1e-2):
    """
    概要:
        Molecule オブジェクトから点群とその対称操作を推定します。
    詳細説明:
        pymatgen.symmetry.analyzer.PointGroupAnalyzer を使用して点群記号と対称操作を特定します。
        toleranceとeigen_toleranceは対称性の検出における許容誤差を制御します。
    引数:
        :param mol: 点群を推定する Molecule オブジェクト。
        :type mol: pymatgen.core.Molecule
        :param tolerance: 座標の許容差 (Å)。
        :type tolerance: float
        :param eigen_tolerance: 固有値縮退の許容差。
        :type eigen_tolerance: float
    戻り値:
        :returns: Schoenflies記号と対称操作のリストのタプル。
        :rtype: tuple
    """
    pga = PointGroupAnalyzer(mol, tolerance=tolerance, eigen_tolerance=eigen_tolerance)
    sch = pga.sch_symbol  # e.g., Td, Oh, T, Th, C3v ...
    ops = pga.get_symmetry_operations()
    return sch, ops

# ====== ラベル付け（各ファミリー） ======

def label_c2v(R: np.ndarray) -> str:
    """
    概要:
        C2v点群の対称操作にラベルを付けます。
    詳細説明:
        単位元、鏡映（σv(xz)またはσv'(yz)）、またはC2回転（C2(z)またはC2）に分類します。
    引数:
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, σv(xz), C2）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    nrm = mirror_normal(R)
    if nrm is not None:
        # x軸・y軸により近い方で σv/σv' を判定
        return "σv(xz)" if abs(nrm[1]) >= abs(nrm[0]) else "σv'(yz)"
    axis, th = rot_axis_angle(R)
    if axis is not None and abs(axis[2]) > 0.999 and np.isclose(th, np.pi, atol=1e-2):
        return "C2(z)"
    if np.isclose(th, np.pi, atol=1e-2): return "C2"
    return "?"

def label_cnv(n: int, R: np.ndarray) -> str:
    """
    概要:
        Cnv点群（C3v, C4v, C6v）の対称操作にラベルを付けます。
    詳細説明:
        単位元、Cn^k(z)回転、または鏡映（σvまたはσd）に分類します。
        nの値に基づいて回転軸と鏡映面を判定します。
    引数:
        :param n: Cn回転軸の次数（3, 4, 6）。
        :type n: int
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, C3^1(z), σv）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    nrm = mirror_normal(R)
    if nrm is not None:
        nxy = _vec_to_unit(np.array([nrm[0], nrm[1], 0.0]))
        v_dirs = [np.array([ 1, 0, 0.0]), np.array([ 0, 1, 0.0]),
                  np.array([-1, 0, 0.0]), np.array([ 0,-1, 0.0])]
        d_dirs = [np.array([ 1, 1, 0.0]), np.array([ 1,-1, 0.0]),
                  np.array([-1, 1, 0.0]), np.array([-1,-1, 0.0])]
        vscore = max(abs(np.dot(nxy, _vec_to_unit(v))) for v in v_dirs)
        dscore = max(abs(np.dot(nxy, _vec_to_unit(d))) for d in d_dirs)
        return "σv" if vscore >= dscore else "σd"
    det = np.linalg.det(R)
    if det > 0:
        axis, th = rot_axis_angle(R)
        if axis is not None and abs(axis[2]) > 0.999:
            k = int(round(th / (2*np.pi/n))) % n
            if n in (4,6) and k == n//2: return "C2(z)"
            if k == 0: return "E"
            return f"C{n}^{k}(z)"
        if np.isclose(th, np.pi, atol=1e-2): return "C2"
    return "?"

def label_dn(n: int, R: np.ndarray) -> str:
    """
    概要:
        Dn点群（D3, D4, D6）の対称操作にラベルを付けます。
    詳細説明:
        単位元、Cn^k(z)回転、またはC2'(xy)回転に分類します。
        z軸周りのCn回転とxy面内のC2回転を識別します。
    引数:
        :param n: Cn回転軸の次数（3, 4, 6）。
        :type n: int
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, C3^1(z), C2'(xy)）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    det = np.linalg.det(R)
    axis, th = rot_axis_angle(R) if det > 0 else (None, None)
    if det > 0 and axis is not None and abs(axis[2]) > 0.999:
        k = int(round(th / (2*np.pi/n))) % n
        if n % 2 == 0 and k == n//2: return "C2(z)"
        if k == 0: return "E"
        return f"C{n}^{k}(z)"
    if det > 0 and axis is not None and abs(axis[2]) < 1e-3 and np.isclose(th, np.pi, atol=1e-2):
        return "C2'(xy)"
    return "?"

def label_t(R: np.ndarray) -> str:
    """
    概要:
        T点群の対称操作にラベルを付けます。
    詳細説明:
        単位元、C3回転（<111>方向）、またはC2回転（<110>または<100>方向）に分類します。
    引数:
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, C3, C2）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    if np.linalg.det(R) < 0: return "?"
    axis, th = rot_axis_angle(R)
    if axis is None: return "?"
    th_deg = np.degrees(th)
    if (np.isclose(th_deg, 120.0, atol=1.0) or np.isclose(th_deg, 240.0, atol=1.0)) and align_score(axis, U111) > 0.99:
        return "C3" if th_deg < 180.0 else "C3^2"
    if np.isclose(th_deg, 180.0, atol=1.0) and (align_score(axis, U110) > 0.99 or align_score(axis, U100) > 0.99):
        return "C2"
    return "?"

def label_th(R: np.ndarray) -> str:
    """
    概要:
        Th点群の対称操作にラベルを付けます。
    詳細説明:
        単位元、反転中心 i、T点群の回転操作、S6回転、またはσh鏡映に分類します。
        行列式が負の操作には逆行列の回転軸と角度を評価します。
    引数:
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, i, C3, S6, σh）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    if np.allclose(R, -np.eye(3), atol=1e-8): return "i"
    det = np.linalg.det(R)
    if det > 0:
        return label_t(R)
    Rt = -R
    axis, th = rot_axis_angle(Rt)
    if axis is None: return "?"
    th_deg = np.degrees(th)
    if align_score(axis, U111) > 0.99:
        if np.isclose(th_deg, 60.0, atol=1.0):  return "S6"
        if np.isclose(th_deg, 300.0, atol=1.0): return "S6^5"
    if np.isclose(th_deg, 180.0, atol=1.0) and (align_score(axis, U110) > 0.99 or align_score(axis, U100) > 0.99):
        return "σh"
    return "?"

def label_td(R: np.ndarray) -> str:
    """
    概要:
        Td点群の対称操作にラベルを付けます。
    詳細説明:
        単位元、C3回転、C2回転、S4回転、またはσd鏡映に分類します。
        各操作の回転軸や法線ベクトルが特定の結晶軸方向（<100>, <110>, <111>）にどれだけ一致するかを評価します。
    引数:
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, C3, S4, σd）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    det = np.linalg.det(R)
    nrm = mirror_normal(R)
    if nrm is not None:
        return "σd" if align_score(nrm, U110) > 0.99 else "mirror"
    if det > 0:
        axis, th = rot_axis_angle(R)
        if axis is None: return "?"
        th_deg = np.degrees(th)
        if (np.isclose(th_deg, 120.0, atol=1.0) or np.isclose(th_deg, 240.0, atol=1.0)) and align_score(axis, U111) > 0.99:
            return "C3" if th_deg < 180.0 else "C3^2"
        if np.isclose(th_deg, 180.0, atol=1.0) and (align_score(axis, U100) > 0.99 or align_score(axis, U110) > 0.99):
            return "C2"
        return "?"
    Rt = -R
    axis, th = rot_axis_angle(Rt)
    if axis is not None and align_score(axis, U100) > 0.99:
        th_deg = np.degrees(th)
        if np.isclose(th_deg, 90.0, atol=1.0):  return "S4"
        if np.isclose(th_deg, 270.0, atol=1.0): return "S4^3"
    return "?"

def label_o(R: np.ndarray) -> str:
    """
    概要:
        O点群の対称操作にラベルを付けます。
    詳細説明:
        単位元、C3回転、C4回転、またはC2回転に分類します。
        回転軸が<100>または<110>または<111>方向であるかを評価します。
    引数:
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, C3, C4, C2(⟨100⟩)）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    if np.linalg.det(R) < 0: return "?"
    axis, th = rot_axis_angle(R)
    if axis is None: return "?"
    th_deg = np.degrees(th)
    if (np.isclose(th_deg, 120.0, atol=1.0) or np.isclose(th_deg, 240.0, atol=1.0)) and align_score(axis, U111) > 0.99:
        return "C3" if th_deg < 180.0 else "C3^2"
    if align_score(axis, U100) > 0.99:
        if np.isclose(th_deg, 90.0, atol=1.0):  return "C4"
        if np.isclose(th_deg, 270.0, atol=1.0): return "C4^3"
        if np.isclose(th_deg, 180.0, atol=1.0): return "C2(⟨100⟩)"
    if np.isclose(th_deg, 180.0, atol=1.0) and align_score(axis, U110) > 0.99:
        return "C2(⟨110⟩)"
    return "?"

def label_oh(R: np.ndarray) -> str:
    """
    概要:
        Oh点群の対称操作にラベルを付けます。
    詳細説明:
        単位元、反転中心 i、O点群の回転操作、S4回転、またはS6回転に分類します。
        回転軸が<100>または<111>方向であるかを評価します。
    引数:
        :param R: 対称操作の回転行列。
        :type R: numpy.ndarray
    戻り値:
        :returns: ラベル文字列（例: E, i, C4, S6）。
        :rtype: str
    """
    R = _orth_preserve_det(R)
    if np.allclose(R, np.eye(3), atol=1e-8): return "E"
    if np.allclose(R, -np.eye(3), atol=1e-8): return "i"
    det = np.linalg.det(R)
    if det > 0:
        return label_o(R)
    Rt = -R
    axis, th = rot_axis_angle(Rt)
    if axis is None: return "?"
    th_deg = np.degrees(th)
    if align_score(axis, U100) > 0.99:
        if np.isclose(th_deg, 90.0, atol=1.0):  return "S4"
        if np.isclose(th_deg, 270.0, atol=1.0): return "S4^3"
    if align_score(axis, U111) > 0.99:
        if np.isclose(th_deg, 60.0, atol=1.0):   return "S6"
        if np.isclose(th_deg, 300.0, atol=1.0):  return "S6^5"
    return "?"

# ====== メイン ======

def main():
    """
    概要:
        コマンドライン引数に基づいてXYZ分子の点群を推定し、対称操作をラベル付けして表示します。
    詳細説明:
        入力XYZファイルから分子を読み込み、オプションに応じて質量中心への平行移動や慣性主軸への整列を行います。
        pymatgen.symmetry.analyzer.PointGroupAnalyzer を使用して点群を推定し、
        tkpointgroup ユーティリティとカスタムのラベリング関数を使って対称操作を分類・表示します。
        点群の推定結果と、選択された点群ファミリーに基づくラベル付け結果が標準出力されます。
        --xyz で入力ファイル、--tolerance で座標許容差、--eigen-tol で固有値縮退の許容差を指定します。
        --center で質量中心へ平行移動、--align で慣性主軸整列を行います。
        --dump-ops を指定すると直交化前の行列も表示され、--write-aligned で整列後のXYZファイルを保存できます。
        --assume オプションで分類ファミリーを明示的に指定することも可能です。
    戻り値:
        なし
    """
    ap = argparse.ArgumentParser(
        description="XYZ 分子から点群を推定し、対称操作を Cnv/Dn/T/Td/Th/O/Oh 記号でラベル付け。"
    )
    ap.add_argument("--xyz", required=True, help="入力 XYZ")
    ap.add_argument("--tolerance", type=float, default=0.03, help="座標許容差 [Å]")
    ap.add_argument("--eigen-tol", type=float, default=1e-2, help="固有値縮退の許容（無次元）")
    ap.add_argument("--center", action="store_true", help="質量中心へ平行移動")
    ap.add_argument("--align", action="store_true", help="慣性主軸整列")
    ap.add_argument("--dump-ops", action="store_true", help="直交化前の行列も表示")
    ap.add_argument("--write-aligned", metavar="OUT_XYZ", help="整列後 XYZ を保存")
    ap.add_argument("--assume",
        choices=["C2v","C3v","C4v","C6v","D3","D4","D6","T","Th","Td","O","Oh"],
        help="分類ファミリー（未指定なら推定記号に基づいて自動）")
    args = ap.parse_args()

    mol = Molecule.from_file(args.xyz)

    if args.center or args.align:
        com = center_of_mass(mol)
        R0 = np.eye(3)
        if args.align:
            R0 = inertia_align_matrix(mol)
        mol = apply_rigid_transform(mol, R0, com)

    sch_raw, ops0 = guess_point_group(mol, tolerance=args.tolerance, eigen_tolerance=args.eigen_tolerance)
    sch = pg.normalize_symbol(sch_raw)  # ← 記号を正規化
    # Pymatgen の操作行列 → 直交化・代表スナップ
    ops = [_orth_preserve_det(op.rotation_matrix) for op in ops0]

    # ファミリー自動選択
    assume = args.assume
    if assume is None:
        if sch in {"Td"}:          assume = "Td"
        elif sch in {"Oh"}:        assume = "Oh"
        elif sch in {"O"}:         assume = "O"
        elif sch in {"T"}:         assume = "T"
        elif sch in {"Th"}:        assume = "Th"
        elif "v" in sch:           assume = sch if sch in {"C3v","C4v","C6v","C2v"} else "C2v"
        elif sch.startswith("D"):  assume = sch if sch in {"D3","D4","D6"} else "D3"
        else:
            if "3" in sch:   assume = "C3v"
            elif "4" in sch: assume = "C4v"
            elif "6" in sch: assume = "C6v"
            else:            assume = "C2v"

    # ラベリング
    labeled = []
    for R in ops:
        if assume == "Td":        name = label_td(R)
        elif assume == "O":       name = label_o(R)
        elif assume == "Oh":      name = label_oh(R)
        elif assume == "T":       name = label_t(R)
        elif assume == "Th":      name = label_th(R)
        elif assume == "C2v":     name = label_c2v(R)
        elif assume in ("C3v","C4v","C6v"):
            n = {"C3v":3,"C4v":4,"C6v":6}[assume]
            name = label_cnv(n, R)
        elif assume in ("D3","D4","D6"):
            n = int(assume[1:])
            name = label_dn(n, R)
        else:
            name = "?"
        # フォールバック（任意）：どうしても判別できないときに基本分類を添える
        if name == "?":
            name = pg.classify_label(R)
        labeled.append((name, R))

    # 出力
    print(f"File            : {args.xyz}")
    print(f"Atoms           : {len(mol)}")
    print(f"Guessed PG      : {sch_raw}")
    print(f"Assume family   : {assume}")
    print(f"Group order     : {len(ops)} operations")

    if args.dump_ops:
        print("\n--- Raw rotation matrices (from analyzer) ---")
        for i, op in enumerate(ops0, 1):
            print(f"[{i:02d}]")
            print(np.array2string(op.rotation_matrix, formatter={'float_kind': lambda x: f"{x:8.5f}"}))

    order_key = {
        "E": 0, "i": 0.5,
        # O/Oh
        "C3": 1, "C3^2": 2, "C4": 1.5, "C4^3": 1.6, "C2(⟨100⟩)": 2.5, "C2(⟨110⟩)": 2.6,
        "S4": 3.0, "S4^3": 3.1, "S6": 3.2, "S6^5": 3.3,
        # Td/T/Th 汎用
        "C2": 2.7, "σd": 4.0, "σh": 4.05, "mirror": 4.1,
        # Cnv/Dn
        "C2(z)": 1.2, "C3^1(z)": 1.3, "C3^2(z)": 1.4, "C4^1(z)": 1.5, "C4^2(z)": 1.6,
        "σv": 4.2, "σd(Cnv)": 4.21, "σv(xz)": 4.22, "σv'(yz)": 4.23, "C2'(xy)": 5,
        "?": 9
    }
    def sort_key(item):
        name, _ = item
        return (order_key.get(name, 8), name)

    print("\n--- Canonicalized & labeled operations ---")
    for i, (name, R) in enumerate(sorted(labeled, key=sort_key), 1):
        print(f"[{i:02d}] {name}\n{pretty_matrix(R)}\n")

    if args.write_aligned:
        out = Path(args.write_aligned)
        write_xyz(mol, out, comment=f"aligned; PG≈{sch_raw}; family={assume}")
        print(f"Aligned XYZ written to: {out}")

if __name__ == "__main__":
    main()
    input("\nPress ENTER to terminate>>\n")