"""
Pymatgenライブラリの各種機能を実演するスクリプトです。
本スクリプトは、CIFファイルを読み込み、その構造に対して様々な解析（基本情報、対称性、XRDシミュレーション、密度計算、空間群操作など）を実行します。
Pymatgenの主要なクラスである`Structure`, `Element`, `Composition`, `SpacegroupAnalyzer`, `XRDCalculator`などの利用方法を示します。
関連リンク: :doc:`template_pymatgen_usage`
"""
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行列を見やすいフォーマットで出力します。
Parameters
----------
:param m: 出力する3x3の行列。
:type m: list[list[float]] or numpy.ndarray
"""
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):
"""
構造オブジェクトの共通する詳細情報（格子、対称操作、サイト）を出力するヘルパー関数。
Parameters
----------
:param structure: 詳細情報を表示するPymatgenのStructureオブジェクト。
:type structure: Structure
:param label: 出力セクションのラベル。
:type 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) を計算します。
Parameters
----------
:param site: PymatgenのSiteオブジェクト。
:type site: pymatgen.core.sites.Site
Returns
-------
:returns: 有効占有数の合計とサイトの総質量 (amu) のタプル。
:rtype: Tuple[float, float]
"""
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]:
"""
構造全体の有効原子数、総質量、サイト数を計算します。
Parameters
----------
:param struct: PymatgenのStructureオブジェクト。
:type struct: Structure
Returns
-------
:returns: 有効原子数、総質量 (amu)、サイト数のタプル。
:rtype: 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]:
"""
構造の格子情報、サイト情報、密度をレポートします。
Parameters
----------
:param struct: PymatgenのStructureオブジェクト。
:type struct: Structure
:param label: 出力セクションのラベル。
:type label: str
Returns
-------
:returns: 原子密度 (atoms/ų) と質量密度 (g/cm³) のタプル。
:rtype: 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を受け取り、その対称操作を分類して文字列で返します。
回転行列の行列式と並進ベクトルの有無に基づいて、Identity、Pure Translation、Inversion、
Rotation (or Screw)、Mirror (or Glide) のいずれかに分類します。
Parameters
----------
:param op: 分類するPymatgenのSymmOpオブジェクト。
:type op: SymmOp
:param tol: 浮動小数点比較の許容誤差。
:type tol: float, optional
Returns
-------
:returns: 対称操作の種類を表す文字列。
:rtype: str
"""
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 expand_coordinates(ispg: int, xyz_str: str, rmin: float = 1e-5) -> List[np.ndarray]:
"""
指定された空間群の対称操作で分率座標を展開し、重複を排除して表示します。
Parameters
----------
:param ispg: 空間群の国際数（International Number）。
:type ispg: int
:param xyz_str: 展開する分率座標の文字列 (例: "0.1, 0.2, 0.3")。
:type xyz_str: str
:param rmin: 重複を判断するための最小距離許容値 (トーラス距離)。
:type rmin: float, optional
Returns
-------
:returns: 独立な展開座標のNumPy配列のリスト。
:rtype: List[numpy.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オブジェクトの基本的なプロパティ（格子、空間群、組成）を表示します。
Parameters
----------
:param structure: 解析対象のPymatgenのStructureオブジェクト。
:type 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):
"""
組成に基づいて酸化数と総電荷を推定し、表示します。
Pymatgenの`composition.oxi_state_guesses()`メソッドを使用して、可能な酸化数の組み合わせを試行します。
Parameters
----------
:param structure: 解析対象のPymatgenのStructureオブジェクト。
:type 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):
"""
構造内のサイトを反復処理し、特定の原子種（金属）のサイトを抽出し、その近傍原子を検索します。
Pymatgenの`site.specie.is_metal`プロパティと`structure.get_neighbors()`メソッドを使用します。
Parameters
----------
:param structure: 解析対象のPymatgenのStructureオブジェクト。
:type 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`クラスの操作をテストし、表示します。
`Element`の基本的なプロパティ、`Composition`の還元化学式と`almost_equals`メソッドの利用を示します。
"""
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) パターンを計算し、主要なピークを表示します。
Pymatgenの`XRDCalculator`クラスを使用して、Cu Kα線による回折パターンをシミュレーションします。
Parameters
----------
:param structure: 解析対象のPymatgenのStructureオブジェクト。
:type structure: Structure
"""
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オブジェクトから、元の構造と対称化された構造の詳細情報（格子、対称操作、サイト）を表示します。
`SpacegroupAnalyzer`を使用して構造の対称性を解析し、`_print_structure_details`ヘルパー関数を呼び出します。
Parameters
----------
:param structure: 解析対象のPymatgenのStructureオブジェクト。
:type 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）と密度チェック ---
# --- 解析機能ブロック 8: 空間群情報表示と座標展開のテスト ---
[ドキュメント]
def test_spacegroup_analysis():
"""
Structureオブジェクトに依存しない、特定の空間群 (ITA標準設定) の情報表示と
座標展開のテストを実行します。
`SpaceGroup.from_int_number()`および`SpacegroupAnalyzer`を使用して空間群のプロパティを取得し、
`expand_coordinates()`関数で指定された座標を展開します。
"""
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ファイルを指定でき、デフォルトでは`ZnO.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()