"""
空間群情報表示およびステレオ投影プロットスクリプト
このスクリプトは、pymatgenライブラリを使用して、指定された空間群に関する情報を
表示したり、特定の初期座標に対する等価な位置を展開してステレオ投影図として
プロットしたりする機能を提供します。
モードには以下の4種類があります。
- 'list': 全ての空間群番号とそのHermann-Mauguin記号をリスト表示します。
- 'inf': 指定された空間群の詳細情報（結晶系、点群、全ての対称操作とその分類）を表示します。
- 'stereo': 指定された空間群の対称操作で初期座標を展開し、結果をXY単位格子に
z座標情報を付加したステレオ投影図としてプロットします。
- 'expand': 指定された空間群の対称操作で初期座標を展開し、重複を排除した
一意な座標のリストを表示します。
使用例:
- 全空間群のリスト表示: python space_group.py --mode list
- 空間群194の詳細情報表示: python space_group.py --mode inf --ispg 194
- 空間群194で座標(0.9,0.1,0.15)を展開し、プロット: python space_group.py --mode stereo --ispg 194 --xyz 0.9,0.1,0.15
- 空間群194で座標(0.9,0.1,0.15)を展開し、リスト表示: python space_group.py --mode expand --ispg 194 --xyz 0.9,0.1,0.15
関連リンク:
:doc:`space_group_usage`
"""
import sys
import argparse
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from pymatgen.core.lattice import Lattice
from pymatgen.core.structure import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.symmetry.groups import SpaceGroup
from pymatgen.core.operations import SymmOp
[ドキュメント]
def initialize():
"""
コマンドライン引数を解析するためのArgumentParserを初期化します。
この関数は、スクリプトが受け入れるコマンドライン引数（モード、空間群番号、
初期座標、重複判定の最小距離）を定義し、それらの説明とデフォルト値を設定します。
ヘルプメッセージはRawTextHelpFormatterを使用して整形されます。
:returns: argparse.ArgumentParser: 初期化されたArgumentParserオブジェクト。
"""
parser = argparse.ArgumentParser(
description="Display space group info or plot stereographic projection.",
formatter_class=argparse.RawTextHelpFormatter
)
# mode引数をキーワード引数に変更し、infとlistを追加
parser.add_argument('--mode', choices=['list', 'inf', 'stereo', 'expand'],
help="'list' : List all space group numbers and symbols\n"
"'inf' : Display detailed space group info\n"
"'stereo' : Plot stereographic projection of expanded points\n"
"'expand' : Expand coordinates and display results")
parser.add_argument('--ispg', type=int,
help="Space group number (1 to 230).")
# 初期値を変更
parser.add_argument('--xyz', type=str, default='0.9,0.1,0.15',
help="Coordinates for expansion, e.g., '0.9,0.1,0.15'.")
parser.add_argument('--rmin', type=float, default=1.0e-5,
help="Minimum distance for judging identical sites.")
return parser
[ドキュメント]
def classify_symop_spg(op, tol=1e-5):
"""
pymatgenのSymmOpオブジェクトを受け取り、その対称操作の種類を分類して文字列で返します。
この関数は、SymmOpオブジェクトの回転行列(R)と並進ベクトル(t)を分析し、
恒等操作、純粋な並進、反転、回転、螺旋回転、鏡映、映進などの種類を特定します。
特に回転操作では回転軸のタイプ(x, y, z, または d(diagonal))を、
鏡映/映進では鏡映面のタイプ(mx, my, mz, または d(diagonal))を判別します。
数値比較には`tol`で指定された許容誤差を使用します。
:param op: pymatgen.core.operations.SymmOp: 分類する対称操作オブジェクト。
:param tol: float: 数値比較に使用する許容誤差。デフォルトは1e-5。
:returns: str: 分類された対称操作の種類を示す文字列（例: "Rotation C3 (x-axis)"）。
"""
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):
if is_pure_translation:
return "Identity"
else:
return "Pure Translation"
if np.allclose(R, -np.eye(3), atol=tol):
return "Inversion"
if np.isclose(det, 1.0, atol=tol):
# 回転または螺旋回転
cosθ = (np.trace(R) - 1) / 2
θ = np.arccos(np.clip(cosθ, -1, 1))
# N回回転のNを計算（0除算や微小な角度の場合のNaN/infを避ける）
if θ > tol:
n = int(np.round(2 * np.pi / θ))
else:
n = 1 # 角度がゼロに近い場合はC1とみなす
vals, vecs = np.linalg.eig(R)
# 固有値1に対応する固有ベクトルが回転軸
axis_vec_candidates = vecs[:, np.isclose(vals, 1.0, atol=tol)].real
axis_vec = np.zeros(3)
if axis_vec_candidates.shape[1] > 0:
axis_vec = axis_vec_candidates[:, 0]
if np.linalg.norm(axis_vec) > tol:
axis_vec /= np.linalg.norm(axis_vec)
axis_type = "d" # diagonal
if np.allclose(np.abs(axis_vec), [1,0,0], atol=tol): axis_type = "x"
elif np.allclose(np.abs(axis_vec), [0,1,0], atol=tol): axis_type = "y"
elif np.allclose(np.abs(axis_vec), [0,0,1], atol=tol): axis_type = "z"
if is_pure_translation:
return f"Rotation C{n} ({axis_type}-axis)"
else:
return f"Screw Rotation C{n} ({axis_type}-axis)"
if np.isclose(det, -1.0, atol=tol):
# 鏡映または映進
# R^2 = I ならば鏡映/映進
if np.allclose(R @ R, np.eye(3), atol=tol):
vals, vecs = np.linalg.eig(R)
# 固有値-1に対応する固有ベクトルが鏡映面の法線ベクトル
nrm_candidates = vecs[:, np.isclose(vals, -1.0, atol=tol)].real
nrm = np.zeros(3)
if nrm_candidates.shape[1] > 0:
nrm = nrm_candidates[:, 0]
if np.linalg.norm(nrm) > tol:
nrm /= np.linalg.norm(nrm)
plane_type = "d" # diagonal
if np.allclose(np.abs(nrm), [1,0,0], atol=tol): plane_type = "mx"
elif np.allclose(np.abs(nrm), [0,1,0], atol=tol): plane_type = "my"
elif np.allclose(np.abs(nrm), [0,0,1], atol=tol): plane_type = "mz"
if is_pure_translation:
return f"Mirror (σ) ({plane_type}-plane)"
else:
return f"Glide Reflection (σ) ({plane_type}-plane)"
return "Unknown"
[ドキュメント]
def display_spg_info(ispg=None):
"""
pymatgenを用いて空間群の情報を表示します。
ispgがNoneの場合、全ての空間群番号とそのHermann-Mauguinシンボルをリスト表示します。
ispgが指定された場合、その空間群の詳細情報（Hermann-Mauguinシンボル、結晶系、点群、
および各対称操作の回転行列、並進ベクトル、変換式、分類）を表示します。
:param ispg: int or None: 表示対象の空間群番号 (1-230)。Noneの場合は全空間群のリストを表示します。
:returns: list[tuple[pymatgen.core.operations.SymmOp, str]]: 取得した対称操作オブジェクトとその分類文字列のリスト。
ispgがNoneの場合は空のリストを返します。
"""
if ispg is None:
print("--- All Space Groups (1 to 230) ---")
for i in range(1, 231):
try:
spg_group = SpaceGroup.from_int_number(i)
print(f"[{i:3d}] {spg_group.symbol}")
except Exception as e:
print(f"[{i:3d}] (Could not retrieve info: {e})")
return []
if not (1 <= ispg <= 230):
print("Error: Invalid space group number. Must be between 1 and 230.")
sys.exit(1)
try:
# pymatgenのSpacegroupAnalyzerはStructureオブジェクトを必要とするため、
# ダミーの立方晶系構造を作成して初期化します。
lattice = Lattice.cubic(1)
dummy_structure = Structure.from_spacegroup(
ispg, lattice, ["C"], [[0, 0, 0]]
)
sga = SpacegroupAnalyzer(dummy_structure)
except Exception as e:
print(f"Error: Could not retrieve data for space group {ispg}. {e}")
sys.exit(1)
print(f"--- Detailed Information for Space Group {ispg} ---")
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--- Symmetry Operations ---")
symm_ops = sga.get_symmetry_operations()
print(f"Total operations: {len(symm_ops)}\n")
op_inf = []
for i, op in enumerate(symm_ops):
kind = classify_symop_spg(op)
op_inf.append((op, kind))
print(f"Operation {i + 1}: {kind}")
print(" Rotation Matrix (R):")
matrix_str = "\n".join([f" [{row[0]:7.4f}, {row[1]:7.4f}, {row[2]:7.4f}]" for row in op.rotation_matrix])
print(matrix_str)
print(" Translation Vector (t):")
print(f" [{op.translation_vector[0]:7.4f}, {op.translation_vector[1]:7.4f}, {op.translation_vector[2]:7.4f}]")
transform_eq = format_transform_spg(op.rotation_matrix, op.translation_vector)
print(f" Transformation: {transform_eq}")
print()
return op_inf
[ドキュメント]
def expand_coordinates(ispg, xyz_str, rmin):
"""
指定された空間群(ITA標準設定)の対称操作で分率座標を展開し、重複を排除して表示します。
初期座標`xyz_str`をpymatgenの`SpaceGroup`オブジェクトから取得した全ての対称操作に適用し、
展開された座標を計算します。
各展開座標は[0,1)の範囲に折り返されます。
その後、`rmin`を閾値として、トーラス空間における最短距離を用いて重複する座標を排除し、
一意な座標のリストを生成して表示します。
:param ispg: int: 空間群番号 (1-230)。
:param xyz_str: str: 初期分率座標を表す文字列（例: "0.9,0.1,0.15"）。
:param rmin: float: 重複するサイトと見なすための最小距離。
:returns: list[numpy.ndarray]: 検出された重複しない展開された分率座標のリスト。
"""
# 入力の解釈
try:
initial_point = np.array(list(map(float, xyz_str.split(','))), dtype=float)
if initial_point.shape != (3,):
raise ValueError
except Exception:
print("Error: Invalid format for --xyz. Please use 'x,y,z'.")
sys.exit(1)
# 空間群の対称操作（原点・設定のずれを避けるため SpaceGroup を直接使用）
try:
spg_group = SpaceGroup.from_int_number(ispg)
symm_ops = spg_group.symmetry_ops # List[SymmOp], すべて分率座標系での操作
except Exception as e:
print(f"Error: Could not retrieve symmetry operations for space group {ispg}. {e}")
sys.exit(1)
expanded_points = []
print(f"--- Space Group {ispg} ({spg_group.symbol}), Expanded points from "
f"({initial_point[0]:.4f}, {initial_point[1]:.4f}, {initial_point[2]:.4f}) ---")
# 展開（分率座標 → 分率座標）
for i, op in enumerate(symm_ops):
transformed = op.operate(initial_point) # 分率座標にそのまま適用
# [0,1) に折り返し
transformed = np.mod(transformed, 1.0)
expanded_points.append(transformed)
kind = classify_symop_spg(op)
print(f"Operation {i + 1:02d} ({kind}): "
f"({transformed[0]:7.4f}, {transformed[1]:7.4f}, {transformed[2]:7.4f})")
# 重複排除（トーラス最短距離）
unique_points = []
for new_p in expanded_points:
is_unique = True
for old_p in unique_points:
diff = np.abs(new_p - old_p)
# トーラス空間での最短距離を考慮 (x, y, z が [0,1] の範囲で周期性を持つ)
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"\nTotal expanded points (including duplicates): {len(expanded_points)}")
print(f"Total unique expanded points (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
[ドキュメント]
def plot_unit_cell_projection(ispg, expanded_points):
"""
展開された座標をXY単位格子にプロットし、z値に応じて色とラベルを決定します。
この関数は、与えられた空間群番号と展開された分率座標のリストを受け取り、
これらの点をXY平面上の単位格子内にプロットします。
各点のプロットマーカーの色とラベルは、z座標のパターン（例: z=0, zと1-zのペア、
zとz±0.5のペアなど）に基づいて決定され、結晶学的な等価位置を視覚的に表現します。
マーカーは 'o' 型で、z=0 (または1) の場合は赤、ペアが存在する場合は黄色、
1-zのみが存在する場合は灰色となります。
:param ispg: int: 空間群番号 (1-230)。
:param expanded_points: list[numpy.ndarray]: `expand_coordinates` 関数から得られた展開済み座標のリスト。
:returns: None
"""
if not (1 <= ispg <= 230):
print("Error: space group must be 1–230.")
sys.exit(1)
fig, ax = plt.subplots(figsize=(6,6))
# 単位格子枠
box = np.array([[0,0],[1,0],[1,1],[0,1],[0,0]])
ax.plot(box[:,0], box[:,1], 'k-')
# 座標のグループ化 (XY座標が同じものはZ値をリストにまとめる)
points_by_xy = {}
for p in expanded_points:
x, y, z = p
key = (round(x, 5), round(y, 5)) # 浮動小数点誤差を考慮して丸める
if key not in points_by_xy:
points_by_xy[key] = []
points_by_xy[key].append(z)
# 各xy座標グループをプロット
eps = 1.0e-5 # 浮動小数点比較の許容誤差
i = 0
for (x, y), z_list in points_by_xy.items():
# z_listから重複を除外しソート
z_set = sorted(set(round(z_val, 5) for z_val in z_list))
# 代表値を 1つ選ぶ（例：最小値）。
z0 = float(z_set[0])
print(f"{i+1:03d}: (x,y)=({x:.5f},{y:.5f}) z_set={z_set} :", end="")
has_z_eq_0 = (abs(z0) < eps or abs(z0 - 1.0) < eps) # z=0またはz=1の点か
has_pair = False
has_pair_p05 = False # z+0.5が存在するか
has_pair_m05 = False # z-0.5が存在するか
has_only_1_minus_z = False # 1-zのみが存在するか（z0自身は存在しない）
# 1−z, z±0.5 が同居しているかを z_set 全体で判定
for _z in z_set:
if abs(_z - (1.0 - z0)) < eps and abs(_z - z0) > eps: # _zがz0と異なり、1-z0に等しい
has_pair = True
if abs((_z - z0) - 0.5) < eps:
has_pair_p05 = True
if abs((_z - z0) + 0.5) < eps:
has_pair_m05 = True
# z0が存在しないが1-z0が存在するようなケース
# has_pairがFalseの場合にのみ確認
if not has_pair:
found_z0 = False
for _z in z_set:
if abs(_z - z0) < eps:
found_z0 = True
break
if not found_z0: # z0自体が存在しない場合
for _z in z_set:
if abs(_z - (1.0 - z0)) < eps:
has_only_1_minus_z = True
break
print(f" has_pair={has_pair} has_pair_p05={has_pair_p05} has_pair_m05={has_pair_m05} has_only_1_minus_z={has_only_1_minus_z}")
marker_color = 'white'
label_text = ''
if has_z_eq_0:
marker_color = 'red'
label_text = '0' # Z=0 の点のラベル
elif has_pair: # z と 1-z の両方が存在する
marker_color = 'yellow'
label_text = '+-'
elif has_pair_p05: # z と z+0.5 が存在する
marker_color = 'yellow'
label_text = '+1/2'
elif has_pair_m05: # z と z-0.5 が存在する
marker_color = 'yellow'
label_text = '-1/2'
elif has_only_1_minus_z: # z0 (代表値) はなく、1-z0 のみがあるような場合
marker_color = 'gray'
label_text = '-' # 1-z の点を示す
ax.scatter(x, y, marker='o', s=200, facecolor=marker_color, edgecolor='black')
ax.text(x, y, label_text, ha='center', va='center', color='k', fontsize=12)
i += 1
ax.set_aspect('equal', 'box')
ax.set_xlim(-0.1,1.1)
ax.set_ylim(-0.1,1.1)
ax.set_xlabel('x (fractional)')
ax.set_ylabel('y (fractional)')
spg_group = SpaceGroup.from_int_number(ispg)
ax.set_title(f"Equivalent Positions in Unit Cell (Space Group {ispg} — {spg_group.symbol})")
plt.tight_layout()
plt.show(block=False)
[ドキュメント]
def main():
"""
スクリプトのエントリポイントです。
コマンドライン引数を解析し、指定されたモードに基づいて対応する関数を呼び出します。
モードが指定されていない場合や、`inf`, `stereo`, `expand` モードで`--ispg`が
欠落している場合はエラーメッセージを表示し、スクリプトを終了します。
matplotlibの図が表示されている場合は、ユーザー入力があるまで終了を待ちます。
:returns: None
"""
parser = initialize()
args = parser.parse_args()
if len(sys.argv) == 1:
parser.print_help()
sys.exit(0)
if args.mode == 'list':
display_spg_info(None)
elif args.mode == 'inf':
if args.ispg is None:
print("Error: --ispg is required for 'inf' mode.")
sys.exit(1)
display_spg_info(args.ispg)
elif args.mode == 'stereo':
if args.ispg is None:
print("Error: --ispg is required for 'stereo' mode.")
sys.exit(1)
# expand_coordinatesで座標を取得し、plot_unit_cell_projectionに渡す
expanded_points = expand_coordinates(args.ispg, args.xyz, args.rmin)
plot_unit_cell_projection(args.ispg, expanded_points)
elif args.mode == 'expand':
if args.ispg is None:
print("Error: --ispg is required for 'expand' mode.")
sys.exit(1)
expand_coordinates(args.ispg, args.xyz, args.rmin)
# matplotlibのウィンドウが開いている場合、ユーザーがEnterを押すまで待機する
if 'matplotlib.pyplot' in sys.modules and plt.get_fignums():
input("\nPress ENTER to terminate>>")
if __name__ == '__main__':
main()
# main関数でinput待機をしていない場合も考慮し、スクリプト終了前に再度待機
input("\nPress ENTER to terminate>>\n")