find_point_group.py ダウンロード/コピー
find_point_group.py
find_point_group.py
1# find_point_group.py (refactored to use tkpointgroup.py)
2#!/usr/bin/env python3
3# -*- coding: utf-8 -*-
4
5"""
6概要:
7 XYZ 分子から点群を推定し、対称操作を直交化、整列、ラベル付けして表示します。
8詳細説明:
9 このスクリプトは、XYZ分子ファイルから点群を推定し、検出された対称操作を分類・ラベル付けするツールです。
10 質量中心への平行移動、慣性主軸への整列などの前処理オプションも提供します。
11 pymatgen を使用してXYZファイルを読み込み、PointGroupAnalyzer で点群記号と対称操作を特定します。
12 対称操作の行列はSVDにより直交化され、tkpointgroup ライブラリの関数を用いて代表値にスナップされます。
13 Cnv (n=3,4,6), C2v, Dn (n=3,4,6), T, Td, Th, O, Oh などの点群ファミリーに対応するカスタムラベリング関数により、
14 各対称操作に詳細なSchoenflies記号のラベルが付与されます。
15 Cnv群ではσvとσdの区別、C2v群ではσv(xz)とσv'(yz)の区別、Dn群ではCn^k(z)とC2'(x系)の分類が行われます。
16 必要となるライブラリは pymatgen です。pip install pymatgen でインストールしてください。
17 tkpointgroup.py はこのスクリプトと同じディレクトリ、またはPYTHONPATHに存在する必要があります。
18関連リンク:
19 find_point_group_usage
20"""
21
22# コマンドライン実行例はDocstringではなく、Pythonのコメントとして保持します。
23# python find_point_group.py --xyz h2o.xyz --center --align --dump-ops --write-aligned h2o_aligned.xyz
24# python find_point_group.py --xyz nh3.xyz --center --align --symprec 2e-3 --angle-tol 8
25# python find_point_group.py --xyz ch4.xyz --center --align --dump-ops
26# python find_point_group.py --xyz h2o.xyz --center --align
27# python find_point_group.py --xyz SF6.xyz --center --align
28# python find_point_group.py --xyz C6H6.xyz --center --align
29# 推奨：オプションはそのまま（eigen-tol=1e-2, matrix-tol=1e-1）
30
31import argparse
32import numpy as np
33from pathlib import Path
34from typing import List, Tuple, Optional
35
36from pymatgen.core import Molecule
37from pymatgen.core.periodic_table import Element
38from pymatgen.symmetry.analyzer import PointGroupAnalyzer
39
40import tkpointgroup as pg # ← 追加：共通ライブラリ
41
42# ====== 数値ユーティリティ（pg を用いつつ、このスクリプト特有の処理を保持） ======
43
44def _orth_preserve_det(R: np.ndarray) -> np.ndarray:
45 """
46 概要:
47 行列を直交化し、detの符号を保ちつつ代表的な値にスナップします。
48 詳細説明:
49 SVDを用いて行列Rを直交行列に変換し、元の行列式detの符号を維持するように調整します。
50 最終的には tkpointgroup.snap_matrix を使用して要素を代表値にスナップします。
51 引数:
52 :param R: 直交化する行列。
53 :type R: numpy.ndarray
54 戻り値:
55 :returns: 直交化され、代表値にスナップされた行列。
56 :rtype: numpy.ndarray
57 """
58 U, _, Vt = np.linalg.svd(R)
59 R_ = U @ Vt
60 if np.sign(np.linalg.det(R_)) != np.sign(np.linalg.det(R)):
61 U[:, -1] *= -1
62 R_ = U @ Vt
63 return pg.snap_matrix(R_) # 最終スナップはライブラリを使用
64
65def _vec_to_unit(v: np.ndarray) -> np.ndarray:
66 """
67 概要:
68 ベクトルを単位ベクトルに正規化します。
69 詳細説明:
70 ベクトルのノルムが非常に小さい場合は、ゼロベクトルを返します。
71 引数:
72 :param v: 正規化するベクトル。
73 :type v: numpy.ndarray
74 戻り値:
75 :returns: 単位ベクトル。
76 :rtype: numpy.ndarray
77 """
78 n = np.linalg.norm(v)
79 return v / n if n > 1e-15 else v
80
81def rot_axis_angle(R: np.ndarray) -> Tuple[Optional[np.ndarray], float]:
82 """
83 概要:
84 回転行列 R の回転軸（Noneのときは不定）と回転角(0..π)を返します。
85 詳細説明:
86 回転行列から固有値と固有ベクトルを計算し、固有値1に対応する固有ベクトルを回転軸とします。
87 回転角はトレースから計算します。回転軸が不定の場合はNoneを返します。
88 引数:
89 :param R: 回転行列。
90 :type R: numpy.ndarray
91 戻り値:
92 :returns: 回転軸（Noneの場合あり）と回転角（ラジアン、0からπの範囲）のタプル。
93 :rtype: tuple
94 """
95 R = _orth_preserve_det(R)
96 if np.allclose(R, np.eye(3), atol=1e-10):
97 return None, 0.0
98 vals, vecs = np.linalg.eig(R)
99 idx = np.where(np.isclose(vals.real, 1.0, atol=1e-6))[0]
100 if idx.size == 0:
101 th = float(np.arccos(np.clip((np.trace(R)-1)/2, -1, 1)))
102 return None, th
103 axis = vecs[:, idx[0]].real
104 axis /= (np.linalg.norm(axis) + 1e-15)
105 th = float(np.arccos(np.clip((np.trace(R)-1)/2, -1, 1)))
106 return axis, th
107
108def mirror_normal(R: np.ndarray) -> Optional[np.ndarray]:
109 """
110 概要:
111 鏡映の法線ベクトルを返します（鏡映でない場合は None）。
112 詳細説明:
113 行列式が負でかつ二乗して単位行列となる行列を鏡映と判定します。
114 鏡映の場合、固有値-1に対応する固有ベクトルを法線ベクトルとします。
115 引数:
116 :param R: 検査する行列。
117 :type R: numpy.ndarray
118 戻り値:
119 :returns: 鏡映の法線ベクトル、またはNone。
120 :rtype: Optional[numpy.ndarray]
121 """
122 R = _orth_preserve_det(R)
123 if not (np.linalg.det(R) < 0 and np.allclose(R @ R, np.eye(3), atol=1e-6)):
124 return None
125 vals, vecs = np.linalg.eig(R)
126 idx = np.where(np.isclose(vals.real, -1.0, atol=1e-6))[0]
127 if idx.size == 0:
128 return None
129 n = vecs[:, idx[0]].real
130 n /= (np.linalg.norm(n) + 1e-15)
131 return n
132
133def pretty_matrix(R: np.ndarray) -> str:
134 """
135 概要:
136 行列を整形して文字列として返します。
137 詳細説明:
138 微小な要素を0として丸め、指定された書式で数値文字列を生成します。
139 引数:
140 :param R: 整形する行列。
141 :type R: numpy.ndarray
142 戻り値:
143 :returns: 整形された行列の文字列表現。
144 :rtype: str
145 """
146 Rt = np.where(np.abs(R) < 1e-10, 0.0, R)
147 Rt = np.round(Rt, 6)
148 return np.array2string(Rt, formatter={'float_kind': lambda x: f"{x:9.6f}"})
149
150# 代表方向（規格化済み）
151U100 = [_vec_to_unit(np.array([sx,0,0])) for sx in (1,-1)] + \
152 [_vec_to_unit(np.array([0,sy,0])) for sy in (1,-1)] + \
153 [_vec_to_unit(np.array([0,0,sz])) for sz in (1,-1)]
154
155U110: List[np.ndarray] = []
156for sx in (1,-1):
157 for sy in (1,-1):
158 U110 += [_vec_to_unit(np.array([sx, sy, 0]))]
159for sy in (1,-1):
160 for sz in (1,-1):
161 U110 += [_vec_to_unit(np.array([0, sy, sz]))]
162for sx in (1,-1):
163 for sz in (1,-1):
164 U110 += [_vec_to_unit(np.array([sx, 0, sz]))]
165
166U111 = [_vec_to_unit(np.array([sx, sy, sz]))
167 for sx in (1,-1) for sy in (1,-1) for sz in (1,-1)]
168
169def align_score(axis: np.ndarray, fam: List[np.ndarray]) -> float:
170 """
171 概要:
172 特定の方向群に対する軸の整列スコアを計算します。
173 詳細説明:
174 与えられた軸が方向群のいずれかのベクトルとどれだけ平行に近いかを示すスコアを計算します。
175 最大の内積の絶対値をスコアとします。
176 引数:
177 :param axis: 評価する軸。
178 :type axis: numpy.ndarray
179 :param fam: 比較する方向ベクトル群。
180 :type fam: List[numpy.ndarray]
181 戻り値:
182 :returns: 整列スコア (0.0から1.0)。
183 :rtype: float
184 """
185 axis = _vec_to_unit(axis)
186 return max(abs(float(np.dot(axis, d))) for d in fam)
187
188# ====== 幾何前処理（共通関数はこのスクリプトで保持） ======
189
190def center_of_mass(mol: Molecule) -> np.ndarray:
191 """
192 概要:
193 分子の質量中心を計算します。
194 引数:
195 :param mol: pymatgen.core.Molecule オブジェクト。
196 :type mol: pymatgen.core.Molecule
197 戻り値:
198 :returns: 分子の質量中心座標。
199 :rtype: numpy.ndarray
200 """
201 masses = np.array([float(Element(s.species_string).atomic_mass) for s in mol.sites])
202 coords = np.array([s.coords for s in mol.sites])
203 return (masses[:, None] * coords).sum(axis=0) / masses.sum()
204
205def inertia_align_matrix(mol: Molecule) -> np.ndarray:
206 """
207 概要:
208 分子の慣性主軸に整列するための回転行列を計算します。
209 詳細説明:
210 質量中心を計算し、それに基づいて慣性テンソルを構築します。
211 慣性テンソルの固有ベクトルが慣性主軸となり、これらを基底とする回転行列を生成します。
212 行列式が負になる場合は、主軸の一つを反転して右手系を保ちます。
213 引数:
214 :param mol: pymatgen.core.Molecule オブジェクト。
215 :type mol: pymatgen.core.Molecule
216 戻り値:
217 :returns: 慣性主軸に整列させるための回転行列。
218 :rtype: numpy.ndarray
219 """
220 masses = np.array([float(Element(s.species_string).atomic_mass) for s in mol.sites])
221 coords = np.array([s.coords for s in mol.sites])
222 com = center_of_mass(mol)
223 r = coords - com
224 I = np.zeros((3, 3))
225 for m, v in zip(masses, r):
226 x, y, z = v
227 I += m * np.array([[y*y+z*z, -x*y, -x*z],
228 [-x*y, x*x+z*z, -y*z ],
229 [-x*z, -y*z, x*x+y*y]])
230 _, vecs = np.linalg.eigh(I)
231 R = vecs.copy()
232 if np.linalg.det(R) < 0:
233 R[:, 0] *= -1
234 return R
235
236def apply_rigid_transform(mol: Molecule, R: np.ndarray, t: np.ndarray) -> Molecule:
237 """
238 概要:
239 分子に剛体変換（回転と平行移動）を適用します。
240 詳細説明:
241 分子の原子座標に平行移動 t を適用し、その後に回転行列 R を適用します。
242 引数:
243 :param mol: 変換を適用する pymatgen.core.Molecule オブジェクト。
244 :type mol: pymatgen.core.Molecule
245 :param R: 回転行列。
246 :type R: numpy.ndarray
247 :param t: 平行移動ベクトル。
248 :type t: numpy.ndarray
249 戻り値:
250 :returns: 変換が適用された新しい Molecule オブジェクト。
251 :rtype: pymatgen.core.Molecule
252 """
253 coords = np.array([s.coords for s in mol.sites])
254 new = (coords - t) @ R.T
255 species = [s.species for s in mol.sites]
256 return Molecule(species, new)
257
258def write_xyz(mol: Molecule, path: Path, comment: str = ""):
259 """
260 概要:
261 Molecule オブジェクトをXYZファイルとして書き込みます。
262 引数:
263 :param mol: 書き込む Molecule オブジェクト。
264 :type mol: pymatgen.core.Molecule
265 :param path: 出力ファイルのパス。
266 :type path: pathlib.Path
267 :param comment: XYZファイルの2行目に書き込むコメント。
268 :type comment: str
269 戻り値:
270 なし
271 """
272 with open(path, "w", encoding="utf-8") as f:
273 f.write(f"{len(mol)}\n{comment}\n")
274 for site in mol.sites:
275 el = site.species_string
276 x, y, z = site.coords
277 f.write(f"{el:2s} {x: .6f} {y: .6f} {z: .6f}\n")
278
279# ====== 点群推定 ======
280
281def guess_point_group(mol: Molecule, tolerance: float = 0.03, eigen_tolerance: float = 1e-2):
282 """
283 概要:
284 Molecule オブジェクトから点群とその対称操作を推定します。
285 詳細説明:
286 pymatgen.symmetry.analyzer.PointGroupAnalyzer を使用して点群記号と対称操作を特定します。
287 toleranceとeigen_toleranceは対称性の検出における許容誤差を制御します。
288 引数:
289 :param mol: 点群を推定する Molecule オブジェクト。
290 :type mol: pymatgen.core.Molecule
291 :param tolerance: 座標の許容差 (Å)。
292 :type tolerance: float
293 :param eigen_tolerance: 固有値縮退の許容差。
294 :type eigen_tolerance: float
295 戻り値:
296 :returns: Schoenflies記号と対称操作のリストのタプル。
297 :rtype: tuple
298 """
299 pga = PointGroupAnalyzer(mol, tolerance=tolerance, eigen_tolerance=eigen_tolerance)
300 sch = pga.sch_symbol # e.g., Td, Oh, T, Th, C3v ...
301 ops = pga.get_symmetry_operations()
302 return sch, ops
303
304# ====== ラベル付け（各ファミリー） ======
305
306def label_c2v(R: np.ndarray) -> str:
307 """
308 概要:
309 C2v点群の対称操作にラベルを付けます。
310 詳細説明:
311 単位元、鏡映（σv(xz)またはσv'(yz)）、またはC2回転（C2(z)またはC2）に分類します。
312 引数:
313 :param R: 対称操作の回転行列。
314 :type R: numpy.ndarray
315 戻り値:
316 :returns: ラベル文字列（例: E, σv(xz), C2）。
317 :rtype: str
318 """
319 R = _orth_preserve_det(R)
320 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
321 nrm = mirror_normal(R)
322 if nrm is not None:
323 # x軸・y軸により近い方で σv/σv' を判定
324 return "σv(xz)" if abs(nrm[1]) >= abs(nrm[0]) else "σv'(yz)"
325 axis, th = rot_axis_angle(R)
326 if axis is not None and abs(axis[2]) > 0.999 and np.isclose(th, np.pi, atol=1e-2):
327 return "C2(z)"
328 if np.isclose(th, np.pi, atol=1e-2): return "C2"
329 return "?"
330
331def label_cnv(n: int, R: np.ndarray) -> str:
332 """
333 概要:
334 Cnv点群（C3v, C4v, C6v）の対称操作にラベルを付けます。
335 詳細説明:
336 単位元、Cn^k(z)回転、または鏡映（σvまたはσd）に分類します。
337 nの値に基づいて回転軸と鏡映面を判定します。
338 引数:
339 :param n: Cn回転軸の次数（3, 4, 6）。
340 :type n: int
341 :param R: 対称操作の回転行列。
342 :type R: numpy.ndarray
343 戻り値:
344 :returns: ラベル文字列（例: E, C3^1(z), σv）。
345 :rtype: str
346 """
347 R = _orth_preserve_det(R)
348 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
349 nrm = mirror_normal(R)
350 if nrm is not None:
351 nxy = _vec_to_unit(np.array([nrm[0], nrm[1], 0.0]))
352 v_dirs = [np.array([ 1, 0, 0.0]), np.array([ 0, 1, 0.0]),
353 np.array([-1, 0, 0.0]), np.array([ 0,-1, 0.0])]
354 d_dirs = [np.array([ 1, 1, 0.0]), np.array([ 1,-1, 0.0]),
355 np.array([-1, 1, 0.0]), np.array([-1,-1, 0.0])]
356 vscore = max(abs(np.dot(nxy, _vec_to_unit(v))) for v in v_dirs)
357 dscore = max(abs(np.dot(nxy, _vec_to_unit(d))) for d in d_dirs)
358 return "σv" if vscore >= dscore else "σd"
359 det = np.linalg.det(R)
360 if det > 0:
361 axis, th = rot_axis_angle(R)
362 if axis is not None and abs(axis[2]) > 0.999:
363 k = int(round(th / (2*np.pi/n))) % n
364 if n in (4,6) and k == n//2: return "C2(z)"
365 if k == 0: return "E"
366 return f"C{n}^{k}(z)"
367 if np.isclose(th, np.pi, atol=1e-2): return "C2"
368 return "?"
369
370def label_dn(n: int, R: np.ndarray) -> str:
371 """
372 概要:
373 Dn点群（D3, D4, D6）の対称操作にラベルを付けます。
374 詳細説明:
375 単位元、Cn^k(z)回転、またはC2'(xy)回転に分類します。
376 z軸周りのCn回転とxy面内のC2回転を識別します。
377 引数:
378 :param n: Cn回転軸の次数（3, 4, 6）。
379 :type n: int
380 :param R: 対称操作の回転行列。
381 :type R: numpy.ndarray
382 戻り値:
383 :returns: ラベル文字列（例: E, C3^1(z), C2'(xy)）。
384 :rtype: str
385 """
386 R = _orth_preserve_det(R)
387 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
388 det = np.linalg.det(R)
389 axis, th = rot_axis_angle(R) if det > 0 else (None, None)
390 if det > 0 and axis is not None and abs(axis[2]) > 0.999:
391 k = int(round(th / (2*np.pi/n))) % n
392 if n % 2 == 0 and k == n//2: return "C2(z)"
393 if k == 0: return "E"
394 return f"C{n}^{k}(z)"
395 if det > 0 and axis is not None and abs(axis[2]) < 1e-3 and np.isclose(th, np.pi, atol=1e-2):
396 return "C2'(xy)"
397 return "?"
398
399def label_t(R: np.ndarray) -> str:
400 """
401 概要:
402 T点群の対称操作にラベルを付けます。
403 詳細説明:
404 単位元、C3回転（<111>方向）、またはC2回転（<110>または<100>方向）に分類します。
405 引数:
406 :param R: 対称操作の回転行列。
407 :type R: numpy.ndarray
408 戻り値:
409 :returns: ラベル文字列（例: E, C3, C2）。
410 :rtype: str
411 """
412 R = _orth_preserve_det(R)
413 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
414 if np.linalg.det(R) < 0: return "?"
415 axis, th = rot_axis_angle(R)
416 if axis is None: return "?"
417 th_deg = np.degrees(th)
418 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:
419 return "C3" if th_deg < 180.0 else "C3^2"
420 if np.isclose(th_deg, 180.0, atol=1.0) and (align_score(axis, U110) > 0.99 or align_score(axis, U100) > 0.99):
421 return "C2"
422 return "?"
423
424def label_th(R: np.ndarray) -> str:
425 """
426 概要:
427 Th点群の対称操作にラベルを付けます。
428 詳細説明:
429 単位元、反転中心 i、T点群の回転操作、S6回転、またはσh鏡映に分類します。
430 行列式が負の操作には逆行列の回転軸と角度を評価します。
431 引数:
432 :param R: 対称操作の回転行列。
433 :type R: numpy.ndarray
434 戻り値:
435 :returns: ラベル文字列（例: E, i, C3, S6, σh）。
436 :rtype: str
437 """
438 R = _orth_preserve_det(R)
439 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
440 if np.allclose(R, -np.eye(3), atol=1e-8): return "i"
441 det = np.linalg.det(R)
442 if det > 0:
443 return label_t(R)
444 Rt = -R
445 axis, th = rot_axis_angle(Rt)
446 if axis is None: return "?"
447 th_deg = np.degrees(th)
448 if align_score(axis, U111) > 0.99:
449 if np.isclose(th_deg, 60.0, atol=1.0): return "S6"
450 if np.isclose(th_deg, 300.0, atol=1.0): return "S6^5"
451 if np.isclose(th_deg, 180.0, atol=1.0) and (align_score(axis, U110) > 0.99 or align_score(axis, U100) > 0.99):
452 return "σh"
453 return "?"
454
455def label_td(R: np.ndarray) -> str:
456 """
457 概要:
458 Td点群の対称操作にラベルを付けます。
459 詳細説明:
460 単位元、C3回転、C2回転、S4回転、またはσd鏡映に分類します。
461 各操作の回転軸や法線ベクトルが特定の結晶軸方向（<100>, <110>, <111>）にどれだけ一致するかを評価します。
462 引数:
463 :param R: 対称操作の回転行列。
464 :type R: numpy.ndarray
465 戻り値:
466 :returns: ラベル文字列（例: E, C3, S4, σd）。
467 :rtype: str
468 """
469 R = _orth_preserve_det(R)
470 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
471 det = np.linalg.det(R)
472 nrm = mirror_normal(R)
473 if nrm is not None:
474 return "σd" if align_score(nrm, U110) > 0.99 else "mirror"
475 if det > 0:
476 axis, th = rot_axis_angle(R)
477 if axis is None: return "?"
478 th_deg = np.degrees(th)
479 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:
480 return "C3" if th_deg < 180.0 else "C3^2"
481 if np.isclose(th_deg, 180.0, atol=1.0) and (align_score(axis, U100) > 0.99 or align_score(axis, U110) > 0.99):
482 return "C2"
483 return "?"
484 Rt = -R
485 axis, th = rot_axis_angle(Rt)
486 if axis is not None and align_score(axis, U100) > 0.99:
487 th_deg = np.degrees(th)
488 if np.isclose(th_deg, 90.0, atol=1.0): return "S4"
489 if np.isclose(th_deg, 270.0, atol=1.0): return "S4^3"
490 return "?"
491
492def label_o(R: np.ndarray) -> str:
493 """
494 概要:
495 O点群の対称操作にラベルを付けます。
496 詳細説明:
497 単位元、C3回転、C4回転、またはC2回転に分類します。
498 回転軸が<100>または<110>または<111>方向であるかを評価します。
499 引数:
500 :param R: 対称操作の回転行列。
501 :type R: numpy.ndarray
502 戻り値:
503 :returns: ラベル文字列（例: E, C3, C4, C2(⟨100⟩)）。
504 :rtype: str
505 """
506 R = _orth_preserve_det(R)
507 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
508 if np.linalg.det(R) < 0: return "?"
509 axis, th = rot_axis_angle(R)
510 if axis is None: return "?"
511 th_deg = np.degrees(th)
512 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:
513 return "C3" if th_deg < 180.0 else "C3^2"
514 if align_score(axis, U100) > 0.99:
515 if np.isclose(th_deg, 90.0, atol=1.0): return "C4"
516 if np.isclose(th_deg, 270.0, atol=1.0): return "C4^3"
517 if np.isclose(th_deg, 180.0, atol=1.0): return "C2(⟨100⟩)"
518 if np.isclose(th_deg, 180.0, atol=1.0) and align_score(axis, U110) > 0.99:
519 return "C2(⟨110⟩)"
520 return "?"
521
522def label_oh(R: np.ndarray) -> str:
523 """
524 概要:
525 Oh点群の対称操作にラベルを付けます。
526 詳細説明:
527 単位元、反転中心 i、O点群の回転操作、S4回転、またはS6回転に分類します。
528 回転軸が<100>または<111>方向であるかを評価します。
529 引数:
530 :param R: 対称操作の回転行列。
531 :type R: numpy.ndarray
532 戻り値:
533 :returns: ラベル文字列（例: E, i, C4, S6）。
534 :rtype: str
535 """
536 R = _orth_preserve_det(R)
537 if np.allclose(R, np.eye(3), atol=1e-8): return "E"
538 if np.allclose(R, -np.eye(3), atol=1e-8): return "i"
539 det = np.linalg.det(R)
540 if det > 0:
541 return label_o(R)
542 Rt = -R
543 axis, th = rot_axis_angle(Rt)
544 if axis is None: return "?"
545 th_deg = np.degrees(th)
546 if align_score(axis, U100) > 0.99:
547 if np.isclose(th_deg, 90.0, atol=1.0): return "S4"
548 if np.isclose(th_deg, 270.0, atol=1.0): return "S4^3"
549 if align_score(axis, U111) > 0.99:
550 if np.isclose(th_deg, 60.0, atol=1.0): return "S6"
551 if np.isclose(th_deg, 300.0, atol=1.0): return "S6^5"
552 return "?"
553
554# ====== メイン ======
555
556def main():
557 """
558 概要:
559 コマンドライン引数に基づいてXYZ分子の点群を推定し、対称操作をラベル付けして表示します。
560 詳細説明:
561 入力XYZファイルから分子を読み込み、オプションに応じて質量中心への平行移動や慣性主軸への整列を行います。
562 pymatgen.symmetry.analyzer.PointGroupAnalyzer を使用して点群を推定し、
563 tkpointgroup ユーティリティとカスタムのラベリング関数を使って対称操作を分類・表示します。
564 点群の推定結果と、選択された点群ファミリーに基づくラベル付け結果が標準出力されます。
565 --xyz で入力ファイル、--tolerance で座標許容差、--eigen-tol で固有値縮退の許容差を指定します。
566 --center で質量中心へ平行移動、--align で慣性主軸整列を行います。
567 --dump-ops を指定すると直交化前の行列も表示され、--write-aligned で整列後のXYZファイルを保存できます。
568 --assume オプションで分類ファミリーを明示的に指定することも可能です。
569 戻り値:
570 なし
571 """
572 ap = argparse.ArgumentParser(
573 description="XYZ 分子から点群を推定し、対称操作を Cnv/Dn/T/Td/Th/O/Oh 記号でラベル付け。"
574 )
575 ap.add_argument("--xyz", required=True, help="入力 XYZ")
576 ap.add_argument("--tolerance", type=float, default=0.03, help="座標許容差 [Å]")
577 ap.add_argument("--eigen-tol", type=float, default=1e-2, help="固有値縮退の許容（無次元）")
578 ap.add_argument("--center", action="store_true", help="質量中心へ平行移動")
579 ap.add_argument("--align", action="store_true", help="慣性主軸整列")
580 ap.add_argument("--dump-ops", action="store_true", help="直交化前の行列も表示")
581 ap.add_argument("--write-aligned", metavar="OUT_XYZ", help="整列後 XYZ を保存")
582 ap.add_argument("--assume",
583 choices=["C2v","C3v","C4v","C6v","D3","D4","D6","T","Th","Td","O","Oh"],
584 help="分類ファミリー（未指定なら推定記号に基づいて自動）")
585 args = ap.parse_args()
586
587 mol = Molecule.from_file(args.xyz)
588
589 if args.center or args.align:
590 com = center_of_mass(mol)
591 R0 = np.eye(3)
592 if args.align:
593 R0 = inertia_align_matrix(mol)
594 mol = apply_rigid_transform(mol, R0, com)
595
596 sch_raw, ops0 = guess_point_group(mol, tolerance=args.tolerance, eigen_tolerance=args.eigen_tolerance)
597 sch = pg.normalize_symbol(sch_raw) # ← 記号を正規化
598 # Pymatgen の操作行列 → 直交化・代表スナップ
599 ops = [_orth_preserve_det(op.rotation_matrix) for op in ops0]
600
601 # ファミリー自動選択
602 assume = args.assume
603 if assume is None:
604 if sch in {"Td"}: assume = "Td"
605 elif sch in {"Oh"}: assume = "Oh"
606 elif sch in {"O"}: assume = "O"
607 elif sch in {"T"}: assume = "T"
608 elif sch in {"Th"}: assume = "Th"
609 elif "v" in sch: assume = sch if sch in {"C3v","C4v","C6v","C2v"} else "C2v"
610 elif sch.startswith("D"): assume = sch if sch in {"D3","D4","D6"} else "D3"
611 else:
612 if "3" in sch: assume = "C3v"
613 elif "4" in sch: assume = "C4v"
614 elif "6" in sch: assume = "C6v"
615 else: assume = "C2v"
616
617 # ラベリング
618 labeled = []
619 for R in ops:
620 if assume == "Td": name = label_td(R)
621 elif assume == "O": name = label_o(R)
622 elif assume == "Oh": name = label_oh(R)
623 elif assume == "T": name = label_t(R)
624 elif assume == "Th": name = label_th(R)
625 elif assume == "C2v": name = label_c2v(R)
626 elif assume in ("C3v","C4v","C6v"):
627 n = {"C3v":3,"C4v":4,"C6v":6}[assume]
628 name = label_cnv(n, R)
629 elif assume in ("D3","D4","D6"):
630 n = int(assume[1:])
631 name = label_dn(n, R)
632 else:
633 name = "?"
634 # フォールバック（任意）：どうしても判別できないときに基本分類を添える
635 if name == "?":
636 name = pg.classify_label(R)
637 labeled.append((name, R))
638
639 # 出力
640 print(f"File : {args.xyz}")
641 print(f"Atoms : {len(mol)}")
642 print(f"Guessed PG : {sch_raw}")
643 print(f"Assume family : {assume}")
644 print(f"Group order : {len(ops)} operations")
645
646 if args.dump_ops:
647 print("\n--- Raw rotation matrices (from analyzer) ---")
648 for i, op in enumerate(ops0, 1):
649 print(f"[{i:02d}]")
650 print(np.array2string(op.rotation_matrix, formatter={'float_kind': lambda x: f"{x:8.5f}"}))
651
652 order_key = {
653 "E": 0, "i": 0.5,
654 # O/Oh
655 "C3": 1, "C3^2": 2, "C4": 1.5, "C4^3": 1.6, "C2(⟨100⟩)": 2.5, "C2(⟨110⟩)": 2.6,
656 "S4": 3.0, "S4^3": 3.1, "S6": 3.2, "S6^5": 3.3,
657 # Td/T/Th 汎用
658 "C2": 2.7, "σd": 4.0, "σh": 4.05, "mirror": 4.1,
659 # Cnv/Dn
660 "C2(z)": 1.2, "C3^1(z)": 1.3, "C3^2(z)": 1.4, "C4^1(z)": 1.5, "C4^2(z)": 1.6,
661 "σv": 4.2, "σd(Cnv)": 4.21, "σv(xz)": 4.22, "σv'(yz)": 4.23, "C2'(xy)": 5,
662 "?": 9
663 }
664 def sort_key(item):
665 name, _ = item
666 return (order_key.get(name, 8), name)
667
668 print("\n--- Canonicalized & labeled operations ---")
669 for i, (name, R) in enumerate(sorted(labeled, key=sort_key), 1):
670 print(f"[{i:02d}] {name}\n{pretty_matrix(R)}\n")
671
672 if args.write_aligned:
673 out = Path(args.write_aligned)
674 write_xyz(mol, out, comment=f"aligned; PG≈{sch_raw}; family={assume}")
675 print(f"Aligned XYZ written to: {out}")
676
677if __name__ == "__main__":
678 main()
679 input("\nPress ENTER to terminate>>\n")