analyze_symmetry.py ダウンロード/コピー

analyze_symmetry.py をダウンロード

analyze_symmetry.py
analyze_symmetry.py
  1"""
  2配位子構造の点群対称性を解析し、対称適合線形結合 (SALC) を構築するスクリプト。
  3
  4詳細説明:
  5本スクリプトは、JSONファイルから配位子の構造データを読み込み、点群対称性を自動検出します。
  6検出された点群に基づいて、配位子の対称性を既約表現に分解し、
  7中心金属のd軌道とのハイブリダイゼーションの可否を評価します。
  8さらに、各既約表現に対応する対称適合線形結合 (SALC) を抽出し、その結果を出力します。
  9
 10関連リンク:
 11:doc:`salc_calculation_overview`
 12"""
 13import json
 14import argparse
 15import numpy as np
 16
 17import tkpg
 18from tkpg import core
 19
 20
 21def load_structure_json(path: str) -> dict:
 22    """
 23    指定されたJSONファイルから構造設定を読み込む。
 24
 25    詳細説明:
 26    ファイルパスを引数として受け取り、配位子の位置、d軌道タイプ、解析モード、
 27    許容誤差などの設定を辞書形式で抽出して返します。
 28    欠損しているキーにはデフォルト値が適用されます。
 29
 30    :param path: str: 構造設定を含むJSONファイルのパス。
 31    :returns: dict: 読み込まれた設定値を含む辞書。
 32        主なキーとして`ligands_pos` (numpy.ndarray), `d_orbital` (str), `mode` (str),
 33        `tol_match_geom` (float), `tol_match_D` (float), `salc_eig_tol` (float),
 34        `coeff_tol` (float), `print_coeffs` (bool), `print_salc_thr` (float),
 35        `max_salc_per_irrep` (int | None), `autodetect` (dict) などが含まれます。
 36    """
 37    with open(path, "r", encoding="utf-8") as f:
 38        obj = json.load(f)
 39
 40    ligands_pos = np.array(obj["ligands_pos"], float)
 41    d_orbital = obj.get("d_orbital", "d_xy")
 42    mode = obj.get("mode", "full")
 43
 44    tol = obj.get("tolerances", {})
 45    tol_match_geom = float(tol.get("tol_match_geom", 5e-3))
 46    tol_match_D = float(tol.get("tol_match_D", 5e-3))
 47
 48    salc_cfg = obj.get("salc", {})
 49    cfg = {
 50        "ligands_pos": ligands_pos,
 51        "d_orbital": d_orbital,
 52        "mode": mode,
 53        "tol_match_geom": tol_match_geom,
 54        "tol_match_D": tol_match_D,
 55        "salc_eig_tol": float(salc_cfg.get("salc_eig_tol", 1e-6)),
 56        "coeff_tol": float(salc_cfg.get("coeff_tol", 1e-6)),
 57        "print_coeffs": bool(salc_cfg.get("print_coeffs", True)),
 58        "print_salc_thr": float(salc_cfg.get("print_salc_thr", 1e-3)),
 59        "max_salc_per_irrep": salc_cfg.get("max_salc_per_irrep", None),
 60        "autodetect": obj.get("autodetect", {}),
 61    }
 62    return cfg
 63
 64
 65def center_positions(ligands_pos: np.ndarray) -> np.ndarray:
 66    """
 67    配位子の位置を重心が原点に来るように移動する。
 68
 69    詳細説明:
 70    与えられた配位子の三次元座標群の重心を計算し、各座標から重心を引くことで、
 71    重心を原点に移動した新しい座標群を生成します。
 72
 73    :param ligands_pos: numpy.ndarray: shapeが(N, 3)の配位子座標配列。
 74    :returns: numpy.ndarray: 重心が原点に移動された配位子座標配列。
 75    """
 76    return ligands_pos - np.mean(ligands_pos, axis=0, keepdims=True)
 77
 78
 79def choose_point_group(
 80    plugins: list[dict],
 81    ligands_pos: np.ndarray,
 82    tol_match_geom: float,
 83    autodetect_cfg: dict | None
 84) -> tuple[dict, list[dict]]:
 85    """
 86    複数の点群プラグインを評価し、最も適切な点群を自動選択する。
 87
 88    詳細説明:
 89    各プラグインが提供する`align_guess`と`symmetry_hit_rate`の結果を評価し、
 90    以下のルールに基づいて最適な点群を選択します:
 91    1. `min_rate_strict`設定で定義された閾値`strict_ok`を満たすプラグインを優先します。
 92    2. 残りの候補の中から、ヒット数 (`hits`) が最大であるものを選択します。
 93    3. 次に、群の位数 (`Gorder`) が最大であるものを選択します。
 94    4. 最後に、ヒット率 (`rate`) が最大であるものを選択します。
 95
 96    :param plugins: list[dict]: 点群プラグインのリスト。
 97        各プラグインは辞書形式で、`name` (str), `align_guess` (callable), `build_group` (callable),
 98        `d_irreps` (dict), `irreps` (callable), `irrep_dim` (callable), `irrep_char` (callable) などの情報を含む。
 99    :param ligands_pos: numpy.ndarray: shapeが(N, 3)の配位子座標配列。
100    :param tol_match_geom: float: 幾何学的な一致を判定するための許容誤差。
101    :param autodetect_cfg: dict | None: 自動検出設定を含む辞書。
102        `min_rate_strict` (dict) などのキーを含むことができる。
103    :returns: tuple[dict, list[dict]]:
104        - `dict`: 選択された最適な点群プラグインと評価結果を含む辞書。
105          キーとして`plugin` (dict), `ok_align` (bool), `R_align` (numpy.ndarray),
106          `pos_aligned` (numpy.ndarray), `hits` (int), `rate` (float),
107          `Gorder` (int), `strict_thr` (float), `strict_ok` (bool) を含む。
108        - `list[dict]`: 評価された全ての点群プラグインの結果リスト。
109          各要素は上記の`dict`と同じ構造を持つ。
110    """
111    autodetect_cfg = autodetect_cfg or {}
112    min_rate = autodetect_cfg.get("min_rate_strict", {})  # 例: {"Oh":0.70, "C4v":0.80}
113
114    scored = []
115    for p in plugins:
116        ok, R_align, pos_aligned = p["align_guess"](ligands_pos)
117        group = p["build_group"]()
118        Gorder = len(group)
119
120        if ok:
121#            hits, rate = core.symmetry_hits(group, pos_aligned, tol_match=tol_match_geom)
122            hits, rate = core.symmetry_hit_rate(group, pos_aligned, tol_match=tol_match_geom)
123        else:
124            hits, rate = 0, 0.0
125
126        # strict acceptance flag (rate based: backward compatible)
127        thr = float(min_rate.get(p["name"], 0.0))
128        strict_ok = (rate >= thr)
129
130        scored.append({
131            "plugin": p,
132            "ok_align": ok,
133            "R_align": R_align,
134            "pos_aligned": pos_aligned,
135            "hits": hits,
136            "rate": rate,
137            "Gorder": Gorder,
138            "strict_thr": thr,
139            "strict_ok": strict_ok,
140        })
141
142    # ---- selection rule ----
143    # 1) Prefer strict_ok True (if any)
144    strict = [s for s in scored if s["strict_ok"]]
145    pool = strict if strict else scored
146
147    # 2) Sort by hits desc, then |G| desc, then rate desc
148    pool.sort(key=lambda s: (s["hits"], s["Gorder"], s["rate"]), reverse=True)
149
150    best = pool[0]
151    return best, scored
152
153
154def main() -> None:
155    """
156    コマンドライン引数からJSONファイルを読み込み、配位子構造の点群対称性解析を実行する。
157
158    詳細説明:
159    1. コマンドライン引数（`--infile`）から構造設定JSONファイルのパスを読み込みます。
160       デフォルトは`structure.json`です。
161    2. JSONファイルから設定を読み込み、配位子の位置を重心が原点となるように正規化します。
162    3. `tkpg`プラグインを使用して最適な点群を自動検出します。
163    4. 検出された点群と設定されたd軌道に基づいて、配位子の既約表現への分解を行います。
164    5. 中心金属d軌道とのハイブリダイゼーションの対称性による可否を判定します。
165    6. 対称適合線形結合 (SALC) を抽出し、その係数を出力します。
166
167    :returns: None
168    """
169    ap = argparse.ArgumentParser()
170    ap.add_argument("--infile", default="structure.json",
171                    help="structure JSON file (default: structure.json)")
172    args = ap.parse_args()
173
174    cfg = load_structure_json(args.infile)
175
176    ligands_pos = center_positions(cfg["ligands_pos"])
177    d_orbital = cfg["d_orbital"]
178    mode = cfg["mode"]
179
180    tol_match_geom = cfg["tol_match_geom"]
181    tol_match_D = cfg["tol_match_D"]
182
183    plugins = tkpg.load_plugins()
184    if not plugins:
185        raise RuntimeError("No tkpg plugins found. Ensure tkpg/Oh.py, tkpg/C4v.py exist.")
186
187    best, scored = choose_point_group(
188        plugins, ligands_pos, tol_match_geom, cfg["autodetect"]
189    )
190
191    plugin = best["plugin"]
192    rate = best["rate"]
193    hits = best["hits"]
194    Gorder = best["Gorder"]
195    R_align = best["R_align"]
196    pos_aligned = best["pos_aligned"]
197
198    print(f"\n[Auto-detected point group]  {plugin['name']}   (hits={hits}/{Gorder}, hit_rate={rate:.3f})")
199
200    print("[Candidates]")
201    # show all candidates, sorted by the SAME selection key (but keep strict_ok shown)
202    scored_sorted = sorted(scored, key=lambda s: (s["strict_ok"], s["hits"], s["Gorder"], s["rate"]), reverse=True)
203    for s in scored_sorted:
204        p = s["plugin"]
205        flag = "STRICT" if s["strict_ok"] else "     "
206        thr = s["strict_thr"]
207        print(f"  {p['name']:>4s} : hits={s['hits']:>2d}/{s['Gorder']:<2d}  rate={s['rate']:.3f}  "
208              f"ok_align={s['ok_align']}  {flag} (thr={thr:.2f})")
209
210    # d_zx alias
211    d_orbital_use = "d_xz" if d_orbital == "d_zx" else d_orbital
212    dmap = plugin["d_irreps"]
213    if d_orbital_use not in dmap:
214        raise ValueError(f"d orbital '{d_orbital}' not defined for point group {plugin['name']}")
215
216    d_irrep = dmap[d_orbital_use]
217    print("\n[Central metal d orbital]")
218    print(f"  Orbital : {d_orbital}")
219    print(f"  Irrep   : {d_irrep}")
220
221    # Build basis
222    frames = core.compute_local_frames(pos_aligned)
223    basis  = core.build_ligand_basis(len(pos_aligned), frames, mode=mode)
224
225    group = plugin["build_group"]()
226    irreps = plugin["irreps"]()
227    irrep_dim = plugin["irrep_dim"]
228    irrep_char = plugin["irrep_char"]
229
230    # Decompose
231    Gamma = core.build_reducible_characters(group, pos_aligned, frames, basis, mode=mode, tol_match=tol_match_D)
232    coeffs = core.decompose_generic(Gamma, group, irreps, irrep_char)
233
234    if cfg["print_coeffs"]:
235        print(f"\n[Ligand reducible rep decomposition ({plugin['name']})]")
236        for ir in irreps:
237            if coeffs.get(ir, 0.0) > cfg["coeff_tol"]:
238                print(f"  {ir}: {coeffs[ir]:.3f}")
239
240    # Hybridization check
241    print("\n[Hybridization check]")
242    if coeffs.get(d_irrep, 0.0) > cfg["coeff_tol"]:
243        print(f"  ✔ Hybridization allowed by symmetry: ligand contains {d_irrep}")
244    else:
245        print(f"  ✘ Symmetry forbidden: ligand does NOT contain {d_irrep}")
246
247    # SALCs: show ALL irreps
248    print(f"\n[All ligand SALCs ({plugin['name']})]  mode={mode}")
249    any_printed = False
250    for ir in irreps:
251        mult = coeffs.get(ir, 0.0)
252        if mult <= cfg["coeff_tol"]:
253            continue
254
255        salcs = core.extract_salc(
256            ir, group,
257            irrep_dim, irrep_char,
258            pos_aligned, frames, basis, mode=mode,
259            tol_match=tol_match_D, tol=cfg["salc_eig_tol"]
260        )
261
262        if len(salcs) == 0:
263            print(f"\n[SALCs for {ir}]  (expected mult≈{mult:.3f}, but none extracted; try lowering salc_eig_tol)")
264            continue
265
266        any_printed = True
267        core.print_salc(
268            ir, salcs, basis,
269            thr=cfg["print_salc_thr"],
270            max_salc=cfg["max_salc_per_irrep"]
271        )
272
273    if not any_printed:
274        print("  (No SALCs extracted. Try increasing tol_match_D or lowering salc_eig_tol.)")
275
276
277if __name__ == "__main__":
278    main()