#!/usr/bin/env python
# infer_symmetry_dof.py
"""
概要:
CIFファイルから結晶構造の格子および原子座標の自由度を推測します。
詳細説明:
pymatgenライブラリを使用して、与えられたCIFファイルから結晶構造を読み込みます。
このスクリプトは、空間群、結晶系、および各独立な原子サイトの情報を分析し、
格子パラメータの自由度と制約、ならびに各原子サイトの分数座標における
独立な自由度を判定します。
原子座標の自由度は、代表サイトの座標を微小に摂動させ、
構造の対称性が維持されるかを検査することによって推測されます。
主な機能:
- 結晶系に基づく格子パラメータの自由度と制約の決定。
- 独立な原子サイトの同定と、そのWyckoff記号、多重度の表示。
- 各独立原子サイトの分数座標における固定値による制約の推測。
- 各独立原子サイトの分数座標における摂動テストによる自由度の推測。
- 空間群操作による等価座標の展開表示。
"""
import argparse
import numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.symmetry.groups import SpaceGroup
[ドキュメント]
def wrap01(x):
"""
概要:
座標値を0から1の範囲に正規化します。
詳細説明:
入力された数値または配列 x の各要素を、1を法とする剰余演算により
0.0以上1.0未満の範囲に変換します。
引数:
:param x: 正規化する座標値または座標配列。
:type x: numpy.ndarray or list or float
戻り値:
:returns: 0以上1未満に正規化された座標値の配列。
:rtype: numpy.ndarray
"""
return np.array(x, dtype=float) % 1.0
[ドキュメント]
def frac_dist(a, b):
"""
概要:
周期境界条件下の2つの分数座標間の距離を計算します。
詳細説明:
2つの分数座標 a と b の差を計算し、その結果を0.5を引いてから1.0で割った剰余を取り、
再度0.5を引くことで、周期境界条件における最短距離差を求めます。
その後、その差のノルム（ユークリッド距離）を返します。
引数:
:param a: 1つ目の分数座標。
:type a: numpy.ndarray or list or float
:param b: 2つ目の分数座標。
:type b: numpy.ndarray or list or float
戻り値:
:returns: 周期境界条件下のユークリッド距離。
:rtype: float
"""
a = np.array(a, dtype=float)
b = np.array(b, dtype=float)
d = (a - b + 0.5) % 1.0 - 0.5
return np.linalg.norm(d)
[ドキュメント]
def close_frac(a, b, tol=1e-5):
"""
概要:
周期境界条件下で2つの分数座標が十分に近似しているかを判定します。
詳細説明:
frac_dist 関数を使用して、2つの分数座標 a と b の間の周期的な距離を計算します。
この距離が与えられた許容誤差 tol 未満であれば、座標は近似していると判断されます。
引数:
:param a: 1つ目の分数座標。
:type a: numpy.ndarray or list or float
:param b: 2つ目の分数座標。
:type b: numpy.ndarray or list or float
:param tol: 許容誤差。
:type tol: float
戻り値:
:returns: 近似していればTrue、そうでなければFalse。
:rtype: bool
"""
return frac_dist(a, b) < tol
[ドキュメント]
def close_mod1(a, b, tol=1e-5):
"""
概要:
1を法として2つの数値が十分に近似しているかを判定します。
詳細説明:
2つの数値 a と b の差を計算し、周期境界条件（1を法とする）を考慮した上で、
その絶対値が許容誤差 tol 未満であるかを判定します。
引数:
:param a: 1つ目の数値。
:type a: float
:param b: 2つ目の数値。
:type b: float
:param tol: 許容誤差。
:type tol: float
戻り値:
:returns: 近似していればTrue、そうでなければFalse。
:rtype: bool
"""
d = abs((a - b + 0.5) % 1.0 - 0.5)
return d < tol
[ドキュメント]
def is_fixed_value(v, values=(0.0, 0.25, 1/3, 0.5, 2/3, 0.75), tol=1e-5):
"""
概要:
指定された値が、既知の固定値のいずれかに近いかを判定します。
詳細説明:
入力された値 v が、あらかじめ定義された固定値のリスト values のいずれかに、
許容誤差 tol の範囲内で1を法として近いかを調べます。
引数:
:param v: 判定する値。
:type v: float
:param values: 比較対象の固定値のタプル。
:type values: tuple
:param tol: 許容誤差。
:type tol: float
戻り値:
:returns: (bool, float or None) 固定値に近ければTrueとその値、そうでなければFalseとNone。
:rtype: tuple
"""
for val in values:
if close_mod1(v, val, tol):
return True, val
return False, None
[ドキュメント]
def unique_frac_coords(coords, tol=1e-8):
"""
概要:
周期境界条件下で重複する分数座標を除去し、ユニークな座標のリストを返します。
詳細説明:
入力された分数座標のリスト coords を走査し、各座標が既にユニークなリストに
含まれているかどうかを close_frac 関数で判定します。
含まれていなければリストに追加し、最終的にユニークな座標のみのリストを返します。
引数:
:param coords: 分数座標のリスト。
:type coords: list of numpy.ndarray
:param tol: 座標の重複を判定するための許容誤差。
:type tol: float
戻り値:
:returns: ユニークな分数座標のリスト。
:rtype: list of numpy.ndarray
"""
unique = []
for c in coords:
c = wrap01(c)
if not any(close_frac(c, u, tol) for u in unique):
unique.append(c)
return unique
[ドキュメント]
def expand_site_by_spacegroup(site_frac, sg_symbol, tol=1e-8):
"""
概要:
指定された空間群記号に基づいて、単一の分数座標から等価なサイトの全座標を展開します。
詳細説明:
pymatgen.symmetry.groups.SpaceGroup クラスを使用して、指定された空間群記号に対応する
対称操作のリストを取得します。これらの対称操作を基準となる分数座標 site_frac に適用し、
生成されたすべての等価サイトの座標を wrap01 で正規化します。
最後に unique_frac_coords 関数を使用して重複する座標を除去し、ユニークなサイトのリストを返します。
引数:
:param site_frac: 基準となる分数座標。
:type site_frac: numpy.ndarray
:param sg_symbol: 空間群記号。
:type sg_symbol: str
:param tol: 重複座標の判定に使用する許容誤差。
:type tol: float
戻り値:
:returns: 空間群によって展開されたユニークな分数座標のリスト。
:rtype: list of numpy.ndarray
"""
sg = SpaceGroup(sg_symbol)
coords = []
for op in sg.symmetry_ops:
coords.append(wrap01(op.operate(site_frac)))
return unique_frac_coords(coords, tol=tol)
[ドキュメント]
def get_lattice_dof_from_crystal_system(crystal_system):
"""
概要:
結晶系に基づいて格子自由度（自由なパラメータと制約）を返します。
詳細説明:
入力された結晶系名（例: triclinic, monoclinic）に対応する、
独立な格子パラメータのリスト（free）と、それらの間の制約のリスト（constraints）を辞書形式で返します。
引数:
:param crystal_system: 結晶系名。
:type crystal_system: str
戻り値:
:returns: 格子自由度に関する情報（"free" と "constraints" のキーを持つ辞書）。
:rtype: dict
例外:
:raises ValueError: 未知の結晶系が指定された場合。
"""
cs = crystal_system.lower()
if cs == "triclinic":
return {
"free": ["a", "b", "c", "alpha", "beta", "gamma"],
"constraints": [],
}
if cs == "monoclinic":
return {
"free": ["a", "b", "c", "beta"],
"constraints": ["alpha = 90", "gamma = 90"],
}
if cs == "orthorhombic":
return {
"free": ["a", "b", "c"],
"constraints": ["alpha = beta = gamma = 90"],
}
if cs == "tetragonal":
return {
"free": ["a", "c"],
"constraints": ["a = b", "alpha = beta = gamma = 90"],
}
if cs in ("trigonal", "hexagonal"):
return {
"free": ["a", "c"],
"constraints": ["a = b", "alpha = beta = 90", "gamma = 120"],
}
if cs == "cubic":
return {
"free": ["a"],
"constraints": ["a = b = c", "alpha = beta = gamma = 90"],
}
raise ValueError(f"Unknown crystal system: {crystal_system}")
[ドキュメント]
def get_symmetry_signature(structure, symprec, angle_tolerance):
"""
概要:
結晶構造から対称性情報を抽出し、対称性シグネチャを生成します。
詳細説明:
pymatgen.symmetry.analyzer.SpacegroupAnalyzer を使用して、入力された結晶構造の
空間群記号、空間群番号、結晶系を判定します。
また、対称化された構造から等価サイト群を抽出し、各サイト群の種、Wyckoff記号、多重度、
および代表サイトの分数座標を収集して、サイト情報のリストを作成します。
これらの情報をまとめて対称性シグネチャとして辞書形式で返します。
引数:
:param structure: 入力結晶構造。
:type structure: pymatgen.core.Structure
:param symprec: 対称性検出の距離許容誤差。
:type symprec: float
:param angle_tolerance: 対称性検出の角度許容誤差。
:type angle_tolerance: float
戻り値:
:returns: 空間群記号、空間群番号、結晶系、サイト情報を含む対称性シグネチャ。
:rtype: dict
"""
sga = SpacegroupAnalyzer(
structure,
symprec=symprec,
angle_tolerance=angle_tolerance,
)
symm = sga.get_symmetrized_structure()
site_info = []
for i, sites in enumerate(symm.equivalent_sites):
site0 = sites[0]
site_info.append({
"species": site0.species_string,
"wyckoff": symm.wyckoff_symbols[i],
"multiplicity": len(sites),
"frac": wrap01(site0.frac_coords),
})
return {
"sg_symbol": sga.get_space_group_symbol(),
"sg_number": sga.get_space_group_number(),
"crystal_system": sga.get_crystal_system(),
"site_info": site_info,
}
[ドキュメント]
def same_symmetry_signature(ref, test):
"""
概要:
2つの対称性シグネチャが同じであるかを判定します。
詳細説明:
2つの対称性シグネチャ ref と test を比較し、
空間群番号、独立サイト群の数、および各サイト群の種、Wyckoff記号、多重度が
全て一致するかを検証します。比較の際には、サイト群の情報をソートして順序に依存しない比較を行います。
引数:
:param ref: 比較対象の基準対称性シグネチャ。
:type ref: dict
:param test: 比較対象のテスト対称性シグネチャ。
:type test: dict
戻り値:
:returns: 両方のシグネチャが同じであればTrue、そうでなければFalse。
:rtype: bool
"""
if ref["sg_number"] != test["sg_number"]:
return False
if len(ref["site_info"]) != len(test["site_info"]):
return False
ref_sig = sorted(
(x["species"], x["wyckoff"], x["multiplicity"])
for x in ref["site_info"]
)
test_sig = sorted(
(x["species"], x["wyckoff"], x["multiplicity"])
for x in test["site_info"]
)
return ref_sig == test_sig
[ドキュメント]
def infer_value_based_constraints(frac, tol=1e-5):
"""
概要:
分数座標から、既知の固定値による制約と、座標成分間の等価性の制約を推測します。
詳細説明:
入力された3次元分数座標 frac の各成分(x, y, z)が、
is_fixed_value 関数で定義された既知の固定値のいずれかに近いかを判定します。
固定値に近い場合は fixed 辞書にその値を記録します。
また、固定値でない座標成分間で、1を法として互いに近いペアがあれば、
それらを equal_pairs リストに追加します。
引数:
:param frac: 3次元分数座標。
:type frac: numpy.ndarray
:param tol: 値の比較に使用する許容誤差。
:type tol: float
戻り値:
:returns: (dict, list) 固定値に固定された座標成分の辞書と、等しいと判定された座標成分のペアのリスト。
:rtype: tuple
"""
labels = ["x", "y", "z"]
frac = wrap01(frac)
fixed = {}
for i, lab in enumerate(labels):
ok, val = is_fixed_value(frac[i], tol=tol)
if ok:
fixed[lab] = val
equal_pairs = []
for i in range(3):
for j in range(i + 1, 3):
if labels[i] not in fixed and labels[j] not in fixed:
if close_mod1(frac[i], frac[j], tol):
equal_pairs.append((labels[i], labels[j]))
return fixed, equal_pairs
[ドキュメント]
def find_site_indices_for_group(structure, group_sites, tol=1e-5):
"""
概要:
pymatgen.SymmetrizedStructure から得られた等価サイト群の各サイトに対応する、
元のStructure内のインデックスを見つけます。
詳細説明:
SymmetrizedStructure から抽出された等価サイト群 group_sites 内の各サイトについて、
元の Structure 内で、同じ化学種を持ち、かつ分数座標が許容誤差 tol の範囲内で
一致するサイトを探します。見つかったサイトのインデックスを記録し、
一度使用したインデックスは再度使用しないように管理します。
引数:
:param structure: 元の結晶構造。
:type structure: pymatgen.core.Structure
:param group_sites: 等価サイト群のリスト。
:type group_sites: list of pymatgen.core.Site
:param tol: サイト位置の比較に使用する許容誤差。
:type tol: float
戻り値:
:returns: 各等価サイトに対応するstructure内のインデックスのリスト。
:rtype: list of int
例外:
:raises RuntimeError: 等価サイトが元の構造内で見つからなかった場合。
"""
indices = []
used = set()
for gs in group_sites:
target_frac = wrap01(gs.frac_coords)
found = None
for i, site in enumerate(structure):
if i in used:
continue
if site.species != gs.species:
continue
if close_frac(wrap01(site.frac_coords), target_frac, tol):
found = i
break
if found is None:
raise RuntimeError(
f"Could not match equivalent site {gs.species_string} "
f"at {target_frac} to original structure."
)
indices.append(found)
used.add(found)
return indices
[ドキュメント]
def match_coords_to_indices(structure, indices, new_coords, tol=1e-4):
"""
概要:
元の等価サイト群 indices と、摂動後に展開した new_coords を対応付けます。
詳細説明:
元の等価サイト群のインデックス indices に対応するサイトの分数座標と、
摂動後に空間群展開された新しい座標 new_coords の間で、
周期境界条件下での距離が最も近いものから順に貪欲にマッチングさせます。
これにより、座標の変化が最小となるように、新しい座標を元のサイト順に並べ替えます。
引数:
:param structure: 元の結晶構造。
:type structure: pymatgen.core.Structure
:param indices: 等価サイト群の元のインデックスのリスト。
:type indices: list of int
:param new_coords: 摂動後に空間群展開された分数座標のリスト。
:type new_coords: list of numpy.ndarray
:param tol: 距離計算に使用する許容誤差。
:type tol: float
戻り値:
:returns: 元のindicesの順序に対応付けられた新しい分数座標のリスト。
:rtype: list of numpy.ndarray
"""
remaining = list(new_coords)
matched = []
for idx in indices:
old_frac = wrap01(structure[idx].frac_coords)
distances = [frac_dist(old_frac, c) for c in remaining]
j = int(np.argmin(distances))
matched.append(wrap01(remaining[j]))
remaining.pop(j)
return matched
[ドキュメント]
def perturb_equivalent_group_coordinate(
structure,
group_indices,
representative_frac,
coord_index,
delta,
sg_symbol,
tol=1e-8,
):
"""
概要:
代表サイトの特定座標を微小変位させ、その結果として生成される等価サイト群の新しい座標を返します。
詳細説明:
指定された代表サイト representative_frac の coord_index で示される座標成分を
微小量 delta だけ変位させます。その後、この変位した代表サイトから
指定された空間群 sg_symbol を用いて等価サイト群全体を再展開します。
もし再展開されたサイトの数が元の等価サイト群 group_indices の数と異なる場合、
Noneを返します。そうでない場合は、 match_coords_to_indices 関数を用いて
新しい座標群を元のサイトの順序に対応付け、
元の構造 structure 内の該当サイトを新しい座標で置き換えた新しい構造を返します。
引数:
:param structure: 元の結晶構造。
:type structure: pymatgen.core.Structure
:param group_indices: 摂動対象の等価サイト群のインデックスリスト。
:type group_indices: list of int
:param representative_frac: 等価サイト群の代表サイトの分数座標。
:type representative_frac: numpy.ndarray
:param coord_index: 摂動する座標のインデックス (0=x, 1=y, 2=z)。
:type coord_index: int
:param delta: 座標を摂動させる微小量。
:type delta: float
:param sg_symbol: 空間群記号。
:type sg_symbol: str
:param tol: 座標比較に使用する許容誤差。
:type tol: float
戻り値:
:returns: 摂動後の新しい構造。サイト数が合わない場合はNone。
:rtype: pymatgen.core.Structure or None
"""
rep2 = wrap01(representative_frac.copy())
rep2[coord_index] += delta
rep2 = wrap01(rep2)
expanded = expand_site_by_spacegroup(rep2, sg_symbol, tol=tol)
if len(expanded) != len(group_indices):
return None
new_coords = match_coords_to_indices(
structure=structure,
indices=group_indices,
new_coords=expanded,
tol=tol,
)
s2 = structure.copy()
for idx, frac in zip(group_indices, new_coords):
s2.replace(
idx,
s2[idx].species,
frac,
coords_are_cartesian=False,
)
return s2
[ドキュメント]
def analyze_independent_sites(structure, symprec, angle_tolerance, tol, delta):
"""
概要:
結晶構造内の独立な原子サイトの自由度を解析します。
詳細説明:
pymatgen.symmetry.analyzer.SpacegroupAnalyzer を使用して構造を対称化し、
独立な原子サイトの情報を取得します。
各等価サイト群について、以下の分析を行います。
1. infer_value_based_constraints 関数を用いて、既知の固定値（0.0, 0.25など）や
座標成分間の等価性に基づく制約を推測します。
2. perturb_equivalent_group_coordinate 関数を用いて、代表サイトの各分数座標
（x, y, z）を微小に変位させ、その際に構造の対称性が維持されるかを確認します。
対称性が維持されればその座標は「自由」と判定され、維持されなければ「固定」と判定されます。
これらの結果をまとめて、各独立サイトの解析結果のリストを返します。
引数:
:param structure: 入力結晶構造。
:type structure: pymatgen.core.Structure
:param symprec: 対称性検出の距離許容誤差。
:type symprec: float
:param angle_tolerance: 対称性検出の角度許容誤差。
:type angle_tolerance: float
:param tol: 座標比較に使用する許容誤差。
:type tol: float
:param delta: 摂動に使用する微小量。
:type delta: float
戻り値:
:returns: (dict, list) 構造の対称性シグネチャと、各独立サイトの解析結果のリスト。
:rtype: tuple
"""
ref = get_symmetry_signature(structure, symprec, angle_tolerance)
sga = SpacegroupAnalyzer(
structure,
symprec=symprec,
angle_tolerance=angle_tolerance,
)
symm = sga.get_symmetrized_structure()
labels = ["x", "y", "z"]
results = []
for group_id, sites in enumerate(symm.equivalent_sites):
site0 = sites[0]
frac = wrap01(site0.frac_coords)
wyckoff = symm.wyckoff_symbols[group_id]
group_indices = find_site_indices_for_group(
structure,
sites,
tol=tol,
)
fixed, equal_pairs = infer_value_based_constraints(frac, tol=tol)
perturb_free = []
perturb_fixed = []
for coord_index, lab in enumerate(labels):
s2 = perturb_equivalent_group_coordinate(
structure=structure,
group_indices=group_indices,
representative_frac=frac,
coord_index=coord_index,
delta=delta,
sg_symbol=ref["sg_symbol"],
tol=tol,
)
if s2 is None:
ok = False
else:
try:
test = get_symmetry_signature(
s2,
symprec=symprec,
angle_tolerance=angle_tolerance,
)
ok = same_symmetry_signature(ref, test)
except Exception:
ok = False
if ok:
perturb_free.append(lab)
else:
perturb_fixed.append(lab)
expanded0 = expand_site_by_spacegroup(
frac,
ref["sg_symbol"],
tol=tol,
)
results.append({
"group_id": group_id,
"site_indices": group_indices,
"species": site0.species_string,
"wyckoff": wyckoff,
"multiplicity": len(sites),
"frac": frac,
"fixed_by_value": fixed,
"equal_by_value": equal_pairs,
"free_by_perturb": perturb_free,
"fixed_by_perturb": perturb_fixed,
"expanded_coords": expanded0,
})
return ref, results
[ドキュメント]
def print_coords(coords, indent=" "):
"""
概要:
分数座標のリストを整形して出力します。
詳細説明:
入力された分数座標のリスト coords の各座標を、指定されたインデント indent を付けて、
小数点以下8桁までの書式で標準出力に出力します。
引数:
:param coords: 出力する分数座標のリスト。
:type coords: list of numpy.ndarray
:param indent: 各行の先頭に追加するインデント文字列。
:type indent: str
戻り値:
:returns: なし。
:rtype: None
"""
for i, c in enumerate(coords):
c = wrap01(c)
print(
f"{indent}op[{i:02d}] : "
f"({c[0]:.8f}, {c[1]:.8f}, {c[2]:.8f})"
)
[ドキュメント]
def main():
"""
概要:
コマンドライン引数を解析し、結晶構造の格子および原子座標の自由度を推測して出力します。
詳細説明:
argparseモジュールを使用してコマンドライン引数（CIFファイルパス、symprec、angle_tolerance、
tol、delta、print_expandedオプション）を解析します。
指定されたCIFファイルからpymatgen.core.Structureオブジェクトを読み込み、
analyze_independent_sites 関数で独立な原子サイトの自由度を分析します。
get_lattice_dof_from_crystal_system 関数で結晶系に基づく格子自由度を取得します。
これらの分析結果を整形して標準出力に表示します。
特に、摂動のdelta値がsymprec値よりも小さい場合に警告メッセージを表示します。
引数:
:param なし:
:type なし:
戻り値:
:returns: なし。
:rtype: None
"""
parser = argparse.ArgumentParser(
description=(
"Infer lattice and atomic coordinate degrees of freedom "
"from CIF using pymatgen."
)
)
parser.add_argument("cif", help="Input CIF file")
parser.add_argument("--symprec", type=float, default=1e-4)
parser.add_argument("--angle-tolerance", type=float, default=5.0)
parser.add_argument("--tol", type=float, default=1e-6)
parser.add_argument("--delta", type=float, default=1e-3)
parser.add_argument(
"--print-expanded",
type=int,
default=1,
choices=[0, 1],
help="Print expanded equivalent coordinates by space group.",
)
args = parser.parse_args()
if args.delta <= args.symprec:
print("WARNING: delta should usually be larger than symprec.")
print(f" symprec = {args.symprec:g}")
print(f" delta = {args.delta:g}")
print(" recommended: delta >= 10 * symprec")
print()
structure = Structure.from_file(args.cif)
ref, site_results = analyze_independent_sites(
structure=structure,
symprec=args.symprec,
angle_tolerance=args.angle_tolerance,
tol=args.tol,
delta=args.delta,
)
lattice_dof = get_lattice_dof_from_crystal_system(ref["crystal_system"])
print("====================================")
print("Symmetry")
print("====================================")
print(f"Space group : {ref['sg_symbol']} ({ref['sg_number']})")
print(f"Crystal system : {ref['crystal_system']}")
print()
print("====================================")
print("Lattice degrees of freedom")
print("====================================")
print("Free lattice parameters:")
print(" " + ", ".join(lattice_dof["free"]))
if lattice_dof["constraints"]:
print("Constraints:")
for c in lattice_dof["constraints"]:
print(f" {c}")
else:
print("Constraints: none")
print()
print("====================================")
print("Independent atomic sites")
print("====================================")
for r in site_results:
frac = r["frac"]
print(f"[{r['group_id']}] species : {r['species']}")
print(f" site indices : {r['site_indices']}")
print(f" Wyckoff : {r['wyckoff']}")
print(f" multiplicity : {r['multiplicity']}")
print(
" frac coords : "
f"({frac[0]:.8f}, {frac[1]:.8f}, {frac[2]:.8f})"
)
if r["fixed_by_value"]:
s = ", ".join(
f"{k}={v:.8g}" for k, v in r["fixed_by_value"].items()
)
else:
s = "none"
print(f" fixed by value : {s}")
if r["equal_by_value"]:
s = ", ".join(f"{a}={b}" for a, b in r["equal_by_value"])
else:
s = "none"
print(f" equal by value : {s}")
if r["free_by_perturb"]:
s = ", ".join(r["free_by_perturb"])
else:
s = "none"
print(f" free by perturbation : {s}")
if r["fixed_by_perturb"]:
s = ", ".join(r["fixed_by_perturb"])
else:
s = "none"
print(f" fixed by perturbation : {s}")
print(f" estimated coordinate dof : {len(r['free_by_perturb'])}")
if args.print_expanded:
print(" expanded coords by space group:")
print_coords(r["expanded_coords"])
print()
if __name__ == "__main__":
main()