vib_irreps.py ダウンロード/コピー
vib_irreps.py
vib_irreps.py
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3
4"""
5点群の指標表を用いて分子の振動既約表現を計算するスクリプト。
6
7概要:
8 指定されたXYZファイルと点群シンボルに基づき、分子の振動モードの既約表現を決定します。
9
10詳細説明:
11 本スクリプトは、分子の対称性を表す点群（Schoenflies表記、例: C2v, D3d, Td, Oh, C3h, C4h, Ch, I, Ih）と
12 分子の構造情報を含むXYZファイルを引数として受け取ります。
13 `pymatgen` が利用可能な環境では、ヘルマン・モーガン表記を内部で使用することがありますが、
14 それ以外の場合や特定の点群（I, Ih, C3hなど）については、既知のクラスサイズと一般位置の仮定に基づく
15 「抽象モード」にフォールバックします。
16
17 計算は、まず全自由度に対応する指標 (Γ_3N) を、固定された原子のカウント（サイトベース）または
18 一般位置の仮定から構築します。その後、分子の並進運動 (Γ_T) と回転運動 (Γ_R) の指標を差し引き、
19 純粋な振動モードに対応する指標 (Γ_v) を導出します。
20 最終的に、このΓ_vを指標の直交性定理を用いて、各既約表現への分解係数（多重度）を決定します。
21
22関連リンク:
23 :doc:`vib_irreps_usage`
24
25使用例:
26 python vib_irreps.py --pg Td --xyz CH4.xyz
27 python vib_irreps.py --pg Oh --xyz SF6.xyz
28 python vib_irreps.py --pg D6h --xyz C6H6.xyz
29 python vib_irreps.py --pg I --xyz C60.xyz
30"""
31
32
33import sys
34import argparse
35import numpy as np
36from typing import Dict, List
37from collections import defaultdict
38
39import tkpointgroup as pg # ★ 共通ライブラリ
40
41# ---------- XYZ loader ----------
42def load_coords_from_xyz(path: str) -> np.ndarray:
43 """
44 指定されたXYZファイルから原子の座標を読み込みます。
45
46 概要:
47 XYZファイルから原子の座標をNumPy配列として読み込みます。
48
49 詳細説明:
50 ファイル内の空行やコメント行、または座標情報を含まない行をスキップし、
51 数値データのみを解析してNumPy配列として返します。
52 最初の数行はヘッダー（原子数、コメント行）であると想定し、
53 適切な行から座標の読み込みを開始します。
54 もし最初の行が整数でない場合（例えばコメント行で始まる場合）は、
55 行の読み込み開始位置を自動的に調整します。
56
57 :param path: str: 読み込むXYZファイルのパス。
58 :returns: numpy.ndarray: 読み込まれた原子座標の配列 (N, 3)。ファイルが空の場合、空の配列を返します。
59 """
60 coords = []
61 with open(path, "r", encoding="utf-8") as f:
62 lines = [ln.strip() for ln in f if ln.strip()]
63 start = 0
64 try:
65 int(lines[0].split()[0]); start = 2
66 except Exception:
67 start = 0
68 for ln in lines[start:]:
69 toks = ln.replace(",", " ").split()
70 if len(toks) < 4:
71 continue
72 x, y, z = map(float, toks[-3:])
73 coords.append([x, y, z])
74 return np.array(coords, dtype=float)
75
76# ---------- Main ----------
77def main():
78 """
79 スクリプトの主処理を実行し、分子の振動既約表現を計算して表示します。
80
81 概要:
82 コマンドライン引数から点群とXYZファイルパスを受け取り、
83 分子の振動モードの既約表現を計算して結果を表示します。
84
85 詳細説明:
86 1. コマンドライン引数（点群シンボル、XYZファイルパス、許容誤差）を解析します。
87 2. 指定された点群シンボルを正規化し、サポートされている点群であることを確認します。
88 3. XYZファイルから原子座標を読み込み、座標が解析できなかった場合はエラーで終了します。
89 4. 点群の特性表と、点群が「抽象モード」（例: I, Ih, C3h など、操作行列が直接定義されていない高対称点群）
90 であるかどうかに応じて、処理を分岐します。
91 5. **通常モードの場合:**
92 - 実際の対称操作行列 (`ops`) を取得します。
93 - 各操作の分類 (`raw_labels`) とクラスへの集約 (`class_map`) を行います。
94 - Γ_3N（全自由度）、Γ_T（並進）、Γ_R（回転）の指標を計算します。
95 - これらの指標を各対称クラスごとに集約し、Γ_v（振動）の指標を導出します。
96 - 最終的に、Γ_vを既約表現に分解します。
97 6. **抽象モードの場合:**
98 - 事前に定義されたクラスサイズ (`pg.ABSTRACT_CLASS_SIZES`) を使用して、
99 Γ_3N、Γ_T、Γ_Rの指標（特にΓ_3Nは'E'クラスで3N、他は0）を計算します。
100 - これらからΓ_v（振動）の指標を導出します。
101 - Γ_vを既約表現に分解します。この際、クラスサイズと群の位数をオーバーライドして使用します。
102 7. 計算された各指標のキャラクターと最終的な振動既約表現の結果を標準出力に整形して表示します。
103
104 :returns: None: 計算結果を標準出力に表示します。
105 """
106 sup = ", ".join(sorted(pg.PG_CHAR_TABLES.keys()))
107 ap = argparse.ArgumentParser(
108 description="Vibrational irreps from XYZ and a point group (Schoenflies-like)."
109 )
110 ap.add_argument("--pg", required=True, help=f"Point group (Schoenflies-like). Supported: {sup}")
111 ap.add_argument("--xyz", required=True, help="XYZ file path")
112 ap.add_argument("--tol", type=float, default=1e-5, help="Tolerance (Å) for fixed atoms")
113 args = ap.parse_args()
114
115 symbol = pg.normalize_symbol(args.pg)
116 if symbol not in pg.PG_CHAR_TABLES:
117 print(f"Error: point group '{args.pg}' not supported. Supported: {sorted(pg.PG_CHAR_TABLES.keys())}")
118 sys.exit(1)
119
120 coords = load_coords_from_xyz(args.xyz)
121 if coords.size == 0:
122 print("Error: no coordinates parsed from XYZ.")
123 sys.exit(1)
124
125 tbl = pg.character_table(symbol)
126 classes = tbl["classes"]
127 abstract = symbol in pg.ABSTRACT_CLASS_SIZES
128
129 if not abstract:
130 # --- 通常モード: 実操作行列を使う ---
131 ops = pg.group_ops(symbol) # 3x3
132 raw_labels = [pg.classify_op_for_table(symbol, R) for R in ops]
133 classes, class_map = pg.class_aggregation(symbol, raw_labels)
134
135 chi_3N_lab = pg.gamma_3N_characters(coords, ops, raw_labels, tol=args.tol)
136 chi_T_lab = pg.gamma_trans_characters(raw_labels, class_map)
137 chi_R_lab = pg.gamma_rot_characters(raw_labels, class_map)
138
139 chi_3N = pg.reduce_to_classes(chi_3N_lab, classes, class_map)
140 chi_T = {c: chi_T_lab.get(c, 0.0) for c in classes}
141 chi_R = {c: chi_R_lab.get(c, 0.0) for c in classes}
142 chi_v = {c: chi_3N[c] - chi_T[c] - chi_R[c] for c in classes}
143
144 mults = pg.decompose_irreps(chi_v, symbol)
145
146 else:
147 # --- 抽象モード: クラスサイズのみから構成（I, Ih, C3h） ---
148 class_sizes = pg.ABSTRACT_CLASS_SIZES[symbol]
149 raw_labels = []
150 for cl, sz in class_sizes.items():
151 raw_labels.extend([cl]*sz)
152 # identity class_map
153 chi_T_lab = pg.gamma_trans_characters(raw_labels, {c:c for c in classes})
154 chi_R_lab = pg.gamma_rot_characters(raw_labels, {c:c for c in classes})
155 chi_T = {c: chi_T_lab.get(c, 0.0) for c in classes}
156 chi_R = {c: chi_R_lab.get(c, 0.0) for c in classes}
157 chi_3N = {c: 0.0 for c in classes}; chi_3N["E"] = 3.0 * coords.shape[0]
158 chi_v = {c: chi_3N[c] - chi_T[c] - chi_R[c] for c in classes}
159
160 mults = pg.decompose_irreps(chi_v, symbol,
161 class_sizes_override=class_sizes,
162 order_override=sum(class_sizes.values()))
163
164 def row(d): return " ".join(f"{cl:>7s}:{d.get(cl,0.0):>6.1f}" for cl in classes)
165
166 print(f"Point group: {symbol} (HM = {pg.schoenflies_to_hm(symbol)}) mode={'abstract' if abstract else 'normal'}")
167 print("Classes :", " ".join(f"{cl:>7s}" for cl in classes))
168 print("Γ_3N chars:", row(chi_3N))
169 print("Γ_T chars:", row(chi_T))
170 print("Γ_R chars:", row(chi_R))
171 print("Γ_v chars:", row(chi_v))
172 print("\nVibrational irreps:", pg.pretty_irreps(symbol, mults))
173
174if __name__ == "__main__":
175 main()
176 input("\nPress ENTER to terminate>>\n")