analyze_symmetry_k.py ダウンロード/コピー
analyze_symmetry_k.py
analyze_symmetry_k.py
1"""
2このモジュールは、結晶構造の対称性を解析し、特に中心金属のd軌道と配位子の対称適応型線形結合 (SALC) のk点における互換性を評価します。
3ユーザーが指定した構造JSONファイルから情報を読み込み、点群の自動検出、SALCの計算、およびk点におけるSALCの振る舞いを解析します。
4
5関連リンク: :doc:`analyze_symmetry_k_usage`
6"""
7
8import json
9import argparse
10import numpy as np
11
12import tkpg
13from tkpg import core
14
15
16def load_structure_json(path):
17 """
18 構造設定をJSONファイルから読み込み、設定辞書を構築します。
19
20 JSONファイルから格子、配位子座標、d軌道情報、解析モード、許容誤差、SALC設定などを読み込みます。
21 `ligands_frac` (分数座標) または `ligands_pos` (デカルト座標) のいずれかが必須です。
22
23 :param path: 設定を読み込むJSONファイルへのパス。
24 :type path: str
25 :returns: 読み込んだ設定情報を含む辞書。
26 :rtype: dict
27 :raises ValueError: JSONファイルに 'ligands_frac' または 'ligands_pos' のいずれも含まれていない場合。
28 """
29 with open(path, "r", encoding="utf-8") as f:
30 obj = json.load(f)
31
32 # --- lattice / positions ---
33 lattice = np.array(obj.get("lattice", np.eye(3)), float)
34
35 if "ligands_frac" in obj:
36 ligands_frac = np.array(obj["ligands_frac"], float)
37 # keep as-is; core will build ligands_pos + periodic metadata
38 ligands_pos = None
39 elif "ligands_pos" in obj:
40 ligands_pos = np.array(obj["ligands_pos"], float)
41 ligands_frac = None
42 else:
43 raise ValueError("structure.json must contain either 'ligands_frac' or 'ligands_pos'.")
44
45 d_orbital = obj.get("d_orbital", "d_xy")
46 mode = obj.get("mode", "full")
47
48 center_ligands = bool(obj.get("center_ligands", False))
49
50 # --- tolerances ---
51 tol = obj.get("tolerances", {})
52 tol_match_geom = float(tol.get("tol_match_geom", 5e-3))
53 tol_match_D = float(tol.get("tol_match_D", 5e-3))
54 tol_frac_site = float(tol.get("tol_frac_site", 1e-6)) # periodic site clustering in frac
55
56 # --- SALC config ---
57 salc_cfg = obj.get("salc", {})
58 cfg = {
59 "lattice": lattice,
60 "ligands_frac": ligands_frac,
61 "ligands_pos": ligands_pos,
62 "center_ligands": center_ligands,
63 "d_orbital": d_orbital,
64 "mode": mode,
65 "tol_match_geom": tol_match_geom,
66 "tol_match_D": tol_match_D,
67 "tol_frac_site": tol_frac_site,
68 "salc_eig_tol": float(salc_cfg.get("salc_eig_tol", 1e-6)),
69 "coeff_tol": float(salc_cfg.get("coeff_tol", 1e-6)),
70 "print_coeffs": bool(salc_cfg.get("print_coeffs", True)),
71 "print_salc_thr": float(salc_cfg.get("print_salc_thr", 1e-3)),
72 "max_salc_per_irrep": salc_cfg.get("max_salc_per_irrep", None),
73 "autodetect": obj.get("autodetect", {}),
74 "kpoints_frac": obj.get("kpoints_frac", [[0.0, 0.0, 0.0]]),
75 "k_tol": float(obj.get("k_tol", 1e-6)),
76 }
77 return cfg
78
79
80def choose_point_group(plugins, ligands_pos, tol_match_geom, autodetect_cfg):
81 """
82 複数の点群プラグインを評価し、配位子の配置に最も適切な点群を自動検出します。
83
84 各プラグインの `align_guess` と `symmetry_hit_rate` を利用して適合度を計算し、
85 ヒット率、ヒット数、群の位数、名前の順でソートして最良のものを選択します。
86 `autodetect_cfg` に `min_rate_strict` が指定されている場合、特定の点群に対して
87 厳密なヒット率の閾値が適用されます。
88
89 :param plugins: 評価する `tkpg` プラグインのリスト。各プラグインは、点群の名前、
90 `align_guess` メソッド、`build_group` メソッドなどを含む辞書である必要があります。
91 :type plugins: list[dict]
92 :param ligands_pos: 配位子の座標 (デカルト座標) のNumPy配列。形状は (N, 3) です。
93 :type ligands_pos: numpy.ndarray
94 :param tol_match_geom: 幾何学的なマッチングに使用する許容誤差。
95 :type tol_match_geom: float
96 :param autodetect_cfg: 自動検出に関する追加設定を含む辞書。
97 `min_rate_strict` (点群名と閾値の辞書) や `tie_eps` (タイブレークの許容誤差)
98 などのキーを含めることができます。`None` の場合、デフォルト値が使用されます。
99 :type autodetect_cfg: dict or None
100 :returns: 以下の2つの要素を含むタプルを返します。
101 - `best`: 検出された最良の点群に関する情報。以下の要素を含むタプルです。
102 - `plugin`: 最良と判断された点群プラグイン辞書。
103 - `rate`: 最良点群のヒット率。
104 - `hits`: 最良点群のヒット数。
105 - `R_align`: 配位子を標準的な軸にアラインメントするための回転行列。
106 - `pos_aligned`: アラインメント後の配位子座標。
107 - `strict_tag`: 厳密な受け入れ基準が適用された場合のタグ文字列。
108 - `scored_pretty`: 評価された全プラグインの情報（整形済み文字列含む）のリスト。
109 各要素は `(rate, hits, ok, p, R, pos, tag)` のタプルです。
110 :rtype: tuple[tuple[dict, float, int, numpy.ndarray, numpy.ndarray, str], list[tuple]]
111 """
112 min_rate = (autodetect_cfg.get("min_rate_strict", {}) if autodetect_cfg else {})
113 eps_tie = float((autodetect_cfg.get("tie_eps", 1e-12)) if autodetect_cfg else 1e-12)
114
115 scored = []
116 for p in plugins:
117 ok, R_align, pos_aligned = p["align_guess"](ligands_pos)
118 group = p["build_group"]()
119 hits, rate = core.symmetry_hit_rate(group, pos_aligned if ok else ligands_pos, tol_match=tol_match_geom)
120 scored.append((rate, hits, ok, p, R_align, pos_aligned))
121
122 scored.sort(key=lambda x: (x[0], x[1], len(x[3]["build_group"]()), x[3]["name"]), reverse=True)
123 best_rate, best_hits, best_ok, best_p, best_R, best_pos = scored[0]
124
125 # strict acceptance if specified
126 thr = float(min_rate.get(best_p["name"], 0.0))
127 strict_tag = ""
128 if best_rate < thr:
129 strict_tag = f" STRICT (thr={thr:.2f})"
130
131 # reformat for printing
132 scored_pretty = []
133 for rate, hits, ok, p, R, pos in scored:
134 thrp = float(min_rate.get(p["name"], 0.0))
135 tag = " OK"
136 if rate < thrp:
137 tag = f" STRICT (thr={thrp:.2f})"
138 scored_pretty.append((rate, hits, ok, p, R, pos, tag))
139
140 return (best_p, best_rate, best_hits, best_R, best_pos, strict_tag), scored_pretty
141
142
143def main():
144 """
145 コマンドライン引数からJSONファイルを読み込み、配位子SALCのk点互換性解析を実行します。
146
147 主な処理フローは以下の通りです。
148 1. 構造JSONを読み込み、解析に必要な設定を抽出します。
149 2. 周期的なメタデータ（`ligands_frac` が指定された場合）を構築し、必要に応じて配位子を中央に配置します。
150 3. `tkpg` プラグインをロードし、`choose_point_group` を使用して最適な点群を自動検出します。
151 4. 検出された点群の標準的な軸の慣例に従って配位子をアラインメントします。
152 5. 中心金属のd軌道に対応する既約表現を特定します。
153 6. 配位子の局所フレームと基底を構築します。
154 7. 点群の既約表現分解（k点独立）を実行し、d軌道とのハイブリダイゼーションの可能性をチェックします。
155 8. 各指定k点について、点群SALCを抽出し、Bloch位相を考慮したk点互換性フィルターを適用し、結果を出力します。
156 """
157 ap = argparse.ArgumentParser()
158 ap.add_argument("--infile", default="structure_k.json", help="structure JSON file (default: structure.json)")
159 args = ap.parse_args()
160
161 cfg = load_structure_json(args.infile)
162
163 # ------------------------------------------------------------
164 # Build geometry + periodic metadata
165 # ------------------------------------------------------------
166 lattice = cfg["lattice"]
167
168 if cfg["ligands_frac"] is not None:
169 meta = core.build_periodic_metadata_from_frac(
170 ligands_frac=cfg["ligands_frac"],
171 lattice=lattice,
172 tol_frac=cfg["tol_frac_site"],
173 )
174 ligands_pos = meta["ligands_pos"] # geometry vectors (cart)
175 site_id = meta["site_id"]
176 cell_T = meta["cell_T"]
177 frac_wrap = meta["frac_wrap"]
178 else:
179 ligands_pos = np.array(cfg["ligands_pos"], float)
180 # non-periodic fallback: every ligand is its own site, T=0
181 site_id = np.arange(len(ligands_pos), dtype=int)
182 cell_T = np.zeros((len(ligands_pos), 3), dtype=int)
183 frac_wrap = None
184
185 # Optional centering (NOT recommended for periodic use)
186 if cfg["center_ligands"]:
187 ligands_pos = core.center_positions(ligands_pos)
188
189 # ------------------------------------------------------------
190 # Load plugins & choose point group
191 # ------------------------------------------------------------
192 plugins = tkpg.load_plugins()
193 if not plugins:
194 raise RuntimeError("No tkpg plugins found. Ensure tkpg/Oh.py, tkpg/C4v.py, etc exist.")
195
196 best, scored = choose_point_group(
197 plugins, ligands_pos, cfg["tol_match_geom"], cfg["autodetect"]
198 )
199 plugin, rate, hits, R_align, pos_aligned, strict_tag = best
200
201 print(f"\n[Auto-detected point group] {plugin['name']} (hits={hits}/{len(plugin['build_group']())}, hit_rate={rate:.3f}){strict_tag}")
202 print("[Candidates]")
203 for r, h, ok, p, *_rest in scored:
204 tag = _rest[-1]
205 print(f" {p['name']:>4s} : hits={h:2d}/{len(p['build_group']()):<2d} rate={r:.3f} ok_align={ok}{tag}")
206
207 # ------------------------------------------------------------
208 # Align geometry for canonical axis conventions
209 # NOTE: periodic metadata (site_id, cell_T) is defined from input frac,
210 # and remains valid; alignment is for symmetry operations only.
211 # ------------------------------------------------------------
212 pos_aligned = (R_align @ ligands_pos.T).T
213
214 # ------------------------------------------------------------
215 # Center d orbital irrep
216 # ------------------------------------------------------------
217 d_orbital = cfg["d_orbital"]
218 d_orbital_use = "d_xz" if d_orbital == "d_zx" else d_orbital
219
220 dmap = plugin["d_irreps"]
221 if d_orbital_use not in dmap:
222 raise ValueError(f"d orbital '{d_orbital}' not defined for point group {plugin['name']}")
223 d_irrep = dmap[d_orbital_use]
224
225 print("\n[Central metal d orbital]")
226 print(f" Orbital : {d_orbital}")
227 print(f" Irrep : {d_irrep}")
228
229 # ------------------------------------------------------------
230 # Build basis
231 # ------------------------------------------------------------
232 mode = cfg["mode"]
233 frames = core.compute_local_frames(pos_aligned)
234 basis = core.build_ligand_basis(len(pos_aligned), frames, mode=mode)
235
236 group = plugin["build_group"]()
237 irreps = plugin["irreps"]()
238 irrep_dim = plugin["irrep_dim"]
239 irrep_char = plugin["irrep_char"]
240
241 # ------------------------------------------------------------
242 # Reducible rep decomposition (point group only; k-independent)
243 # ------------------------------------------------------------
244 Gamma = core.build_reducible_characters(
245 group, pos_aligned, frames, basis, mode=mode, tol_match=cfg["tol_match_D"]
246 )
247 coeffs = core.decompose_generic(Gamma, group, irreps, irrep_char)
248
249 if cfg["print_coeffs"]:
250 print(f"\n[Ligand reducible rep decomposition ({plugin['name']})]")
251 for ir in irreps:
252 if coeffs.get(ir, 0.0) > cfg["coeff_tol"]:
253 print(f" {ir}: {coeffs[ir]:.6f}")
254
255 # Symmetry-only hybridization check
256 print("\n[Hybridization check (symmetry-only)]")
257 if coeffs.get(d_irrep, 0.0) > cfg["coeff_tol"]:
258 print(f" ✔ Allowed by point-group symmetry: ligand contains {d_irrep}")
259 else:
260 print(f" ✘ Forbidden by point-group symmetry: ligand does NOT contain {d_irrep}")
261
262 # ------------------------------------------------------------
263 # SALCs (point-group) -> then k-compatibility filtering per k
264 # ------------------------------------------------------------
265 kpoints = np.array(cfg["kpoints_frac"], float)
266 k_tol = float(cfg["k_tol"])
267
268 print(f"\n[All ligand SALCs with Bloch phase] point_group={plugin['name']} mode={mode}")
269 for ir in irreps:
270 mult = coeffs.get(ir, 0.0)
271 if mult <= cfg["coeff_tol"]:
272 continue
273
274 # Extract point-group SALCs (real vectors)
275 salcs = core.extract_salc(
276 ir, group,
277 irrep_dim, irrep_char,
278 pos_aligned, frames, basis, mode=mode,
279 tol_match=cfg["tol_match_D"],
280 tol=cfg["salc_eig_tol"]
281 )
282
283 if len(salcs) == 0:
284 print(f"\n[SALCs for {ir}] (expected mult≈{mult:.3f}, but none extracted; try lowering salc_eig_tol)")
285 continue
286
287 # For each k: filter by k-compatibility, then print with phase
288 for k_frac in kpoints:
289 k_frac = np.array(k_frac, float)
290 compatible, n_ok = core.filter_salcs_by_k(
291 salcs=salcs,
292 basis=basis,
293 site_id=site_id,
294 cell_T=cell_T,
295 k_frac=k_frac,
296 tol=k_tol
297 )
298 print(f"\n[k-compatibility] irrep={ir} k_frac=({k_frac[0]:.3f},{k_frac[1]:.3f},{k_frac[2]:.3f}) compatible={n_ok}/{len(salcs)}")
299
300 if n_ok == 0:
301 print(f" (No k-compatible SALCs for irrep {ir} at k=({k_frac[0]:.3f},{k_frac[1]:.3f},{k_frac[2]:.3f}))")
302 continue
303
304 core.print_salc(
305 ir, compatible, basis,
306 thr=cfg["print_salc_thr"],
307 max_salc=cfg["max_salc_per_irrep"],
308 k_frac=k_frac,
309 cell_T=cell_T
310 )
311
312
313if __name__ == "__main__":
314 main()