import argparse import sys import numpy as np from pprint import pprint from typing import Tuple, List, Optional from math import sqrt, pi, e # Pymatgenの依存関係は、最新の仕様に合わせてサブモジュールから明示的にインポート from pymatgen.core import Structure, Element, Composition from pymatgen.analysis.diffraction.xrd import XRDCalculator from pymatgen.io.cif import CifParser from pymatgen.symmetry.analyzer import SpacegroupAnalyzer from pymatgen.symmetry.groups import SpaceGroup from pymatgen.core.operations import SymmOp # --- ヘルパー関数 --- def print_matrix(m): """3x3行列を見やすいフォーマットで出力します。""" print(f"| {m[0][0]:8.4f} {m[0][1]:8.4f} {m[0][2]:8.4f} |") print(f"| {m[1][0]:8.4f} {m[1][1]:8.4f} {m[1][2]:8.4f} |") print(f"| {m[2][0]:8.4f} {m[2][1]:8.4f} {m[2][2]:8.4f} |") def _print_structure_details(structure: Structure, label: str): """ 構造オブジェクトの共通する詳細情報（格子、対称操作、サイト）を出力するヘルパー関数。 """ print(f"\n--- {label} 詳細情報 ---") # 1. Structure Details structure_dict = structure.as_dict() print(f"Structure keys: {list(structure_dict.keys())}") # 2. Symmetry Analysis analyzer = SpacegroupAnalyzer(structure) spg_name, spg_num = structure.get_space_group_info() symmetry_ops = analyzer.get_symmetry_operations() symmetry_matrix = [op.as_dict()['matrix'] for op in symmetry_ops] nsym = len(symmetry_matrix) print(f"\n空間群: {spg_name} #{spg_num}") print(f"対称操作の数: {nsym}") # 全ての対称操作の出力は冗長なため、最初の2つのみ表示 for i, matrix in enumerate(symmetry_matrix[:2]): print(f" op #{i+1} (Matrix):") print_matrix(matrix) # 3. Lattice Information print("\n格子情報 (Lattice): ") a, b, c, alpha, beta, gamma = structure.lattice.parameters volume = structure.lattice.volume aij = structure.lattice.matrix print(f" 格子定数 (a, b, c): {a:.4f}, {b:.4f}, {c:.4f} Å") print(f" 角度 (α, β, γ): {alpha:.4f}, {beta:.4f}, {gamma:.4f} °") print(f" 体積 (Volume): {volume:.4f} ų") print(" 格子ベクトル行列 (Matrix):") pprint(aij) # 4. Independent Sites (Symmetrized Structureのみがこの情報を持つ) eq_sites = getattr(structure, "equivalent_sites", None) print("\n独立サイト (Independent Sites):") if eq_sites is not None: # 各等価サイトセットの最初のサイト（独立サイト）の情報を表示 for i, sites in enumerate(eq_sites): # sites[0].speciesはCompositionオブジェクト {Element: occupation} species_str = ", ".join([f"{s.name}({occ:.3f})" for s, occ in sites[0].species.items()]) print(f" {i+1}. 種: {species_str}, 分率座標: {sites[0].frac_coords}") else: # Structure.from_fileで作成した構造は通常この属性を持たない print(" 情報なし (SpacegroupAnalyzerで対称化されていないため)") # 5. All sites (最初の5つのみ表示) print("\n全サイト情報 (最初の5つ):") for i, site in enumerate(structure.sites[:5]): composition = site.species species_info = ", ".join([f"{atom.name}({occ:.3f})" for atom, occ in composition.items()]) print(f" {i+1}. 種: {species_info}, 分率座標: {site.frac_coords}") # --- 単位格子変換・密度計算用のヘルパー関数 (機能 7) --- def site_effective_occupancy_and_mass_amu(site) -> Tuple[float, float]: """サイトの有効占有数 ($\Sigma$ occ) と総質量 (amu) を計算します。""" occ = 0.0 mass = 0.0 # site.speciesは {Element: occupation} の辞書 for sp, amt in site.species.items(): occ += float(amt) amu = getattr(sp, "atomic_mass", 0.0) # atomic_massがNoneの場合に対応 mass += float(amt) * float(amu if amu is not None else 0.0) return occ, mass def structure_effective_counts(struct: Structure) -> Tuple[float, float, int]: """構造全体の有効原子数、総質量、サイト数を計算します。""" eff_atoms = 0.0 total_mass = 0.0 for site in struct.sites: occ, mass = site_effective_occupancy_and_mass_amu(site) eff_atoms += occ total_mass += mass return eff_atoms, total_mass, len(struct.sites) def report_structure_density(struct: Structure, label: str) -> Tuple[float, float]: """構造の格子情報、サイト情報、密度をレポートします。""" a, b, c = struct.lattice.abc alpha, beta, gamma = struct.lattice.angles volume = struct.volume eff_atoms, total_mass_amu, n_sites = structure_effective_counts(struct) # 密度の単位変換定数 # 1 amu = 1.66053906660e-24 g # 1 Å^3 = 1.0e-24 cm^3 total_mass_g = total_mass_amu * 1.66053906660e-24 volume_cm3 = volume * 1.0e-24 # 密度計算 rho_atom = eff_atoms / volume if volume > 0 else float("nan") # atoms/Å^3 rho_mass = total_mass_g / volume_cm3 if volume_cm3 > 0 else float("nan") # g/cm^3 print(f"\n[{label} Structure]") print(f" 格子定数: a={a:.6f}, b={b:.6f}, c={c:.6f}") print(f" 角度: $\\alpha={alpha:.2f}$, $\\beta={beta:.2f}$, $\\gamma={gamma:.2f}$") print(f" 体積 (Volume) : {volume:.6f} ų") print(f" サイト数 (Sites) : {n_sites}") print(f" 有効原子数 ($\\Sigma$ occ) : {eff_atoms:.6f}") print(f" 総質量 (Mass) : {total_mass_amu:.6f} amu") print(f" 原子密度 (Atomic $\\rho$): {rho_atom:.6f} atoms/ų") print(f" 質量密度 (Mass $\\rho$): {rho_mass:.6f} g/cm³") return rho_atom, rho_mass # --- 空間群解析・展開用のヘルパー関数 (機能 8) --- def classify_symop_spg(op: SymmOp, tol=1e-5): """SymmOpを受け取り、その対称操作を分類して文字列で返します。""" R = op.rotation_matrix t = op.translation_vector det = np.linalg.det(R) is_pure_translation = np.allclose(t, [0,0,0], atol=tol) if np.allclose(R, np.eye(3), atol=tol): return "Identity" if is_pure_translation else "Pure Translation" if np.allclose(R, -np.eye(3), atol=tol): return "Inversion" if np.isclose(det, 1.0, atol=tol): # 回転操作 (Rotation or Screw) return "Rotation (or Screw)" # 詳細な分類は冗長なため簡略化 if np.isclose(det, -1.0, atol=tol): # 鏡映操作 (Mirror or Glide) return "Mirror (or Glide)" # 詳細な分類は冗長なため簡略化 return "Unknown" def format_transform_spg(matrix, translation_vector): """回転行列と並進ベクトルを座標変換の数式としてフォーマットします。""" labels = ['x', 'y', 'z'] equations = [] for row_idx, row in enumerate(matrix): terms = [] for col_idx, coeff in enumerate(row): if not np.isclose(coeff, 0): term_str = "" if np.isclose(coeff, 1): term_str = labels[col_idx] elif np.isclose(coeff, -1): term_str = f"-{labels[col_idx]}" else: term_str = f"{coeff:.2f}{labels[col_idx]}" if terms and not term_str.startswith('-'): terms.append(f"+{term_str}") else: terms.append(term_str) t_comp = translation_vector[row_idx] if not np.isclose(t_comp, 0): t_str = f"{t_comp:.3f}" if terms and t_comp > 0: terms.append(f"+{t_str}") elif terms and t_comp < 0: terms.append(t_str) else: terms.append(t_str) equations.append("".join(terms) if terms else "0") return f"({equations[0]}, {equations[1]}, {equations[2]})" def expand_coordinates(ispg: int, xyz_str: str, rmin: float = 1e-5) -> List[np.ndarray]: """ 指定された空間群の対称操作で分率座標を展開し、重複を排除して表示します。 """ # 1. 入力の解釈 initial_point = np.array(list(map(float, xyz_str.split(','))), dtype=float) # 2. 空間群の対称操作を取得 try: spg_group = SpaceGroup.from_int_number(ispg) symm_ops = spg_group.symmetry_ops except Exception as e: print(f"エラー: 空間群 {ispg} の対称操作を取得できませんでした: {e}") return [] expanded_points = [] print(f"\n--- 座標展開 (Space Group {ispg} - {spg_group.symbol}) ---") print(f" 初期点: ({initial_point[0]:.4f}, {initial_point[1]:.4f}, {initial_point[2]:.4f})") # 3. 展開（分率座標 $\\rightarrow$ 分率座標） for i, op in enumerate(symm_ops): transformed = op.operate(initial_point) # [0,1) に折り返し transformed = np.mod(transformed, 1.0) expanded_points.append(transformed) # 変換式と結果の表示 (最初の5つのみ) if i < 5: transform_eq = format_transform_spg(op.rotation_matrix, op.translation_vector) print(f" Op {i + 1:02d}: {transform_eq} -> ({transformed[0]:7.4f}, {transformed[1]:7.4f}, {transformed[2]:7.4f})") # 4. 重複排除（トーラス最短距離） unique_points = [] for new_p in expanded_points: is_unique = True for old_p in unique_points: diff = np.abs(new_p - old_p) # トーラス距離 (周期境界条件) diff = np.where(diff > 0.5, 1.0 - diff, diff) if np.linalg.norm(diff) < rmin: is_unique = False break if is_unique: unique_points.append(new_p) print(f"\n総展開点数 (重複含む): {len(expanded_points)}") print(f"総独立点数 (rmin={rmin}): {len(unique_points)}") print("\n--- 独立な展開座標 (Unique Expanded Points) ---") for i, p in enumerate(unique_points): print(f"{i + 1:02d}: ({p[0]:7.4f}, {p[1]:7.4f}, {p[2]:7.4f})") return unique_points # --- 解析機能ブロック 1: 構造の基本情報 --- def analyze_structure_basic(structure: Structure): """ Structureオブジェクトの基本的なプロパティ（格子、空間群、組成）を表示します。 """ print("\n--- 1. 基本構造情報と辞書出力 ---") print(f"型: {type(structure)}") # as_dict()は完全な情報を含むため、辞書のキー数のみ表示 print(f"辞書形式 (as_dict) のキー数: {len(structure.as_dict())}") print(f"化学式: {structure.formula}") print(f"空間群情報: {structure.get_space_group_info()}") print(f"格子定数: {structure.lattice}") # --- 解析機能ブロック 2: 酸化数と電荷 --- def guess_oxidation_states(structure: Structure): """ 組成に基づいて酸化数と総電荷を推定します。 """ print("\n--- 2. 酸化数推定と電荷情報 ---") try: # 可能な酸化数の組み合わせを取得 ox_guesses = structure.composition.oxi_state_guesses() # oxi_state_guessesは可能な組み合わせのリストを返すため、最初の結果を表示 print(f"推定される酸化数 (一例): {ox_guesses[0] if ox_guesses else 'N/A'}") except ValueError as e: # 酸化数推定器が失敗した場合の処理 print(f"酸化数の推定に失敗しました: {e}") print(f"Structureオブジェクトの総電荷: {structure.charge}") # --- 解析機能ブロック 3: サイトと近傍原子の分析 --- def analyze_sites_and_neighbors(structure: Structure): """ サイトを反復処理し、特定の原子種（金属）のサイトを抽出し、近傍原子を検索します。 """ print("\n--- 3. サイト情報と近傍原子検索 ---") # サイトリストの型と、最初のサイト情報を表示 print(f"サイトリストの型: {type(structure.sites)}") if structure.sites: print(f"最初のサイト情報: {structure.sites[0]}") # 金属サイトのフィルタリング metal_sites = [site for site in structure.sites if site.specie.is_metal] print(f"\n全サイト数: {len(structure.sites)}, 検出された金属サイト数: {len(metal_sites)}") # 近傍原子の検索（金属サイトが存在する場合のみ） if metal_sites: # 最初の金属サイトを中心に、半径 5 Å 以内の近傍原子を検索 target_site = metal_sites[0] # cutoff=5.0を指定 neighbors = structure.get_neighbors(target_site, 5.0) print(f"\n--- {target_site.specie.name}サイトの近傍検索 (R=5.0 Å) ---") print(f"近傍原子数: {len(neighbors)}") # 結果の最初の3つのみ表示 for i, neighbor in enumerate(neighbors[:3]): # neighborオブジェクトには距離 (nn_distance) が含まれている print(f" {i+1}. {neighbor.species_string} - 距離: {neighbor.nn_distance:.4f} Å") else: print("\n近傍検索: 金属サイトが見つからなかったためスキップしました。") # --- 解析機能ブロック 4: 元素と組成の操作（Structure非依存） --- def analyze_element_and_composition(): """ Structureオブジェクトに依存しない、ElementとCompositionの操作をテストします。 """ print("\n--- 4. Elementと組成の操作 ---") print("\n[Element: Si 情報]") si = Element('Si') print(f"電子構造: {si.full_electronic_structure}") print(f"金属か: {si.is_metal}") print("\n[Composition: Ti3O5 操作]") ti3o5 = Composition("Ti3O5") # reduced_formulaは安定したプロパティ print(f"還元化学式: {ti3o5.reduced_formula}") print(f"構成元素: {ti3o5.elements}") # 組成の比較（almost_equals）のテスト print("\n[Composition 比較 (almost_equals)]") # reduced_composition プロパティを使用（多くのバージョンでサポート） ti3o5_norm = ti3o5.reduced_composition ti_non_stoich = Composition('Ti2.98O5.10') ti_non_stoich_norm = ti_non_stoich.reduced_composition print(f" 基準組成: {ti3o5_norm}") print(f" 比較組成: {ti_non_stoich_norm}") # 2つの式がある許容範囲内（デフォルト 1e-4）で等しいか判別 is_equal = ti3o5_norm.almost_equals(ti_non_stoich_norm) print(f" almost_equalsによる比較結果: {is_equal}") # --- 解析機能ブロック 5: XRDシミュレーション --- def simulate_xrd(structure: Structure): """ StructureオブジェクトからX線回折 (XRD) パターンを計算し、主要なピークを表示します。 """ print("\n--- 5. XRDシミュレーション ---") # Cu Kα線 (波長 $\\lambda=1.54184$ Å) を使用するXRDCalculatorを初期化 xrd_calc = XRDCalculator(wavelength='CuKa') # XRDパターンを計算 (2$\theta$ = 20° から 80°) xrd_pattern = xrd_calc.get_pattern(structure, two_theta_range=(20, 80)) print(f"XRDピークの総数 (20°～80°): {len(xrd_pattern.x)}") print("\n--- 主要なXRDピーク (最初の5つ) ---") for i in range(min(5, len(xrd_pattern.x))): two_theta = xrd_pattern.x[i] intensity = xrd_pattern.y[i] hkl_planes = xrd_pattern.hkls[i] # hklsは辞書のリスト（最も強い反射の面指数情報を抽出） hkl_data = hkl_planes[0]['hkl'] # 3指数 (h k l) または 4指数 (h k i l) の両方に対応するため、リストを文字列に変換 hkl_str = " ".join(map(str, hkl_data)) print(f" {i+1}. 2$\\theta$: {two_theta:.3f}°, 強度: {intensity:.2f}%, 面指数: ({hkl_str})") # --- 解析機能ブロック 6: CIF/構造情報の詳細分析（SymmetryAnalyzer） --- def analyze_detailed_cif_info(structure: Structure): """ Structureオブジェクトから、元の構造と対称化された構造の詳細情報（格子、対称操作、サイト）を表示します。 """ print("\n" + "="*50) print("--- 6. CIF/構造情報の詳細分析（SymmetryAnalyzer） ---") print("="*50) # 1. 元の構造の解析 _print_structure_details(structure, "元の入力構造") # 2. 対称化された構造の解析 analyzer = SpacegroupAnalyzer(structure) symmetrized_structure = analyzer.get_symmetrized_structure() # 対称化は成功したが、情報が元の構造と大きく変わらない場合がある if symmetrized_structure: _print_structure_details(symmetrized_structure, "対称化された構造") else: print("\n対称化された構造の取得に失敗しました。") # --- 解析機能ブロック 7: 単位格子変換（Primitive Cell）と密度チェック --- def perform_primitive_conversion(structure: Structure, sym_tol: float = 0.01, density_eps: float = 1.0e-4): """ 元の構造を Primitive Standard Structureに変換し、変換前後で密度が一致するか検証します。 """ print("\n" + "="*50) print("--- 7. 単位格子変換（Primitive Cell）と密度チェック ---") print("="*50) # 1. 元の構造のレポート (密度を記録) print("\n[ステップ 1] 変換前の構造情報と密度計算") rho_atom_orig, rho_mass_orig = report_structure_density(structure, "Original") # 2. Primitive Cellへの変換 try: print("\n[ステップ 2] Primitive Standard Structureへの変換 (SpacegroupAnalyzerを使用)") analyzer = SpacegroupAnalyzer(structure, symprec=sym_tol) # get_primitive_standard_structure() は最も厳密なプリミティブ変換を行う s_conv = analyzer.get_primitive_standard_structure() except Exception as e: # 失敗した場合、フォールバックとして get_primitive_structure() を試みる print(f"警告: Standard Structureの取得に失敗。原始格子への変換を試みます: {e}") s_conv = structure.get_primitive_structure() # 3. 変換後の構造のレポート print("\n[ステップ 3] 変換後の構造情報と密度計算") rho_atom_new, rho_mass_new = report_structure_density(s_conv, "Converted (Primitive)") # 4. 密度の一致確認 print("\n[ステップ 4] 密度の一貫性チェック") # 原子密度のチェック is_atom_rho_consistent = np.isclose(rho_atom_orig, rho_atom_new, rtol=density_eps) if is_atom_rho_consistent: print(f"✅ 原子密度は一致しています (許容誤差 $\\epsilon={density_eps}$)") else: print(f"❌ 警告: 原子密度が一致しませんでした。Original: {rho_atom_orig:.6f}, Converted: {rho_atom_new:.6f}") # 質量密度のチェック is_mass_rho_consistent = np.isclose(rho_mass_orig, rho_mass_new, rtol=density_eps) if is_mass_rho_consistent: print(f"✅ 質量密度は一致しています (許容誤差 $\\epsilon={density_eps}$)") else: print(f"❌ 警告: 質量密度が一致しませんでした。Original: {rho_mass_orig:.6f}, Converted: {rho_mass_new:.6f}") print(f"\n[Info] Primitive Structureの体積は、元の体積の {s_conv.volume / structure.volume:.4f} 倍です。") # --- 解析機能ブロック 8: 空間群情報表示と座標展開のテスト --- def test_spacegroup_analysis(): """ Structureオブジェクトに依存しない、特定の空間群 (ITA標準設定) の情報表示と 座標展開のテストを実行します。 """ print("\n" + "="*50) print("--- 8. 空間群情報表示と座標展開のテスト ---") print("="*50) # テストケース 1: 空間群情報の詳細表示 (例: Pm-3m, No. 221) ispg_info = 221 print(f"\n[テスト 1] 空間群 #{ispg_info} の詳細情報") try: # Pymatgenの機能を使用するため、ダミーの Structure を作成する必要がある lattice = Structure.from_spacegroup( ispg_info, [[1, 0, 0], [0, 1, 0], [0, 0, 1]], ["C"], [[0, 0, 0]] ).lattice dummy_structure = Structure.from_spacegroup( ispg_info, lattice, ["C"], [[0, 0, 0]] # Simple cubic cell size 1.0 ) sga = SpacegroupAnalyzer(dummy_structure) print(f"Hermann-Mauguin Symbol: {sga.get_space_group_symbol()}") print(f"Crystal System: {sga.get_crystal_system()}") print(f"Point Group: {sga.get_point_group_symbol()}") print("\n--- 対称操作 (SymmOp) の種類分類 (最初の5つ) ---") symm_ops = sga.get_symmetry_operations() print(f"Total operations: {len(symm_ops)}\n") for i, op in enumerate(symm_ops[:5]): kind = classify_symop_spg(op) transform_eq = format_transform_spg(op.rotation_matrix, op.translation_vector) print(f" Op {i + 1:02d} ({kind}): {transform_eq}") except Exception as e: print(f"エラー: 空間群 {ispg_info} の情報表示中にエラーが発生しました: {e}") # テストケース 2: 座標展開 (例: Fm-3m, No. 225) ispg_expand = 225 xyz_test = '0.1, 0.2, 0.3' rmin_test = 1e-5 print(f"\n[テスト 2] 空間群 #{ispg_expand} での座標展開テスト") # 展開関数の実行 expanded_points = expand_coordinates(ispg_expand, xyz_test, rmin_test) # グラフ表示はスキップ (plt.show()による停止を避けるため) # --- メイン関数 --- def main(): """ CIFファイルを読み込み、すべての解析関数を順次実行します。 """ # デフォルトのCIFファイルを指定 (テスト実行用) default_cif = "ZnO.cif" parser = argparse.ArgumentParser( description="Pymatgenの機能を実行する構造解析プログラム。", formatter_class=argparse.RawTextHelpFormatter ) parser.add_argument( "cif_file", type=str, nargs='?', # オプション引数として扱う default=default_cif, help=f"解析するCIFファイルのパス。\n(デフォルト: {default_cif} を想定)" ) args = parser.parse_args() print(f"Pymatgen解析プログラムを開始します。対象ファイル: {args.cif_file}\n") # Structureオブジェクトの作成 try: # Structure.from_file() を使用して、CIF/POSCARなどのファイルを読み込む structure = Structure.from_file(args.cif_file) print("✅ Structureオブジェクトの作成に成功しました。\n") # 統合された機能を順次実行 analyze_structure_basic(structure) guess_oxidation_states(structure) analyze_sites_and_neighbors(structure) # 構造オブジェクトに依存しない機能 analyze_element_and_composition() simulate_xrd(structure) analyze_detailed_cif_info(structure) perform_primitive_conversion(structure) # 新しい機能ブロック (空間群情報表示と座標展開のテスト) test_spacegroup_analysis() print("\n--- 全ての機能の実行を完了しました。 ---") except FileNotFoundError: print(f"\nエラー: 指定されたファイルが見つかりません: {args.cif_file}", file=sys.stderr) print("💡 ヒント: 実行ファイルと同じディレクトリにCIFファイルが存在するか確認してください。", file=sys.stderr) sys.exit(1) except Exception as e: print(f"\nエラー: Structureオブジェクトの作成または解析中に予期せぬエラーが発生しました: {e}", file=sys.stderr) sys.exit(1) if __name__ == "__main__": main()