#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
pymatgenを用いた格子変換プログラム（部分占有を考慮）。
このプログラムは、pymatgenライブラリを使用して結晶構造の格子変換を実行します。
部分占有サイトも適切に処理し、変換前後の構造の体積、有効原子数、総質量、原子密度、質量密度を報告します。
特に、密度の整合性を検査し、変換によって物理量が変わっていないことを確認します。
主な機能:
- `prim`: SpacegroupAnalyzerを用いた原始セルへの変換。
- `romb`: 六方晶設定のR-格子を菱面体晶の原始セルに変換。
- `hex`: 菱面体晶設定のR-格子を六方晶の慣用セルに変換。
- `orth`: 中心格子（A/B/C/F/I）を適切な変換で原始セルに変換。
- `MATRIX`: ユーザー定義の変換行列 `'(a,b,c)(d,e,f)(g,h,i)'` または `'(a,b,c,tx)(d,e,f,ty)(g,h,i,tz)'` を適用。
行列のエントリは算術式（例: 1/3, sqrt(2)/2）も使用可能。
- 基底変換は `V' = T @ V`、サイト座標は `f' = f @ inv(T) + t` で行い、重複サイトを結合します。
- 出力はセルの体積、有効原子数（占有率の合計）、総質量、原子密度、質量密度を含みます。
- 密度が指定されたeps（デフォルト1e-4）内で一致することを検証します。
関連リンク: :doc:`convert_cell_usage`
"""
import argparse
import sys
import os
import re
import math
import numpy as np
from typing import Tuple, Optional, List
from pymatgen.core.structure import Structure
from pymatgen.core.lattice import Lattice
from pymatgen.io.cif import CifParser
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
# ------------------------- CLI -------------------------
[ドキュメント]
def terminate():
"""
プログラムの実行を終了させ、ユーザーからの入力を待つ。
詳細説明:
標準入力からENTERキーの入力を受け付けるまで待機し、
その後、プログラムを強制終了します。
"""
input("\nPress ENTER to terminate>>")
exit()
[ドキュメント]
def initialize():
"""
コマンドライン引数を解析し、プログラムの設定を初期化する。
詳細説明:
argparseモジュールを使用して、入力ファイル、変換タイプ、
対称性許容誤差などのコマンドライン引数を定義し、解析します。
:returns: argparse.Namespace: 解析されたコマンドライン引数を含むオブジェクト。
"""
p = argparse.ArgumentParser(description="Lattice conversion program using pymatgen.")
p.add_argument("input_file", type=str, help="Input CIF file path.")
p.add_argument("-c", "--conversion", type=str, default="prim",
help=("Conversion spec: prim | rhomb | hex | orth | "
"MATRIX like '(a,b,c)(d,e,f)(g,h,i)' "
"or '(a,b,c,tx)(d,e,f,ty)(g,h,i,tz)'. "
"Each a,b,c,... can be arithmetic (e.g., 1/2, sqrt(2)/2)."))
p.add_argument("-d", "--direction", type=str, default="OriginalToConverted",
choices=["OriginalToConverted", "ConvertedToOriginal"],
help="If a MATRIX is given, reverse by using its inverse when ConvertedToOriginal.")
p.add_argument("-t", "--sym_tol", type=float, default=0.01,
help="Symmetry tolerance for space group analysis.")
p.add_argument("-x", "--xyz_tol", type=float, default=1.0e-4,
help="Tolerance to merge duplicate fractional sites after basis change.")
p.add_argument("--eps", type=float, default=1.0e-4,
help="Relative tolerance for density equality check.")
return p.parse_args()
# ------------------------- eval-based parser -------------------------
# eval() を安全に実行するためのグローバル変数辞書。
# 許可される組み込み関数やモジュール関数のみをリストアップする。
SAFE_GLOBALS = {
"__builtins__": {},
"sqrt": math.sqrt,
"sin": math.sin,
"cos": math.cos,
"tan": math.tan,
"pi": math.pi,
"e": math.e,
}
[ドキュメント]
def print_matrix(message: str, T: np.ndarray):
"""
与えられた行列を整形して標準出力に表示する。
:param message: str: 行列の前に表示するメッセージ。
:param T: np.ndarray: 表示する行列。
"""
print(message)
for row in T:
print(f" |{row[0]:10.6f} {row[1]:10.6f} {row[2]:10.6f}|" )
[ドキュメント]
def parse_number(expr: str) -> float:
"""
文字列の数値式を安全に評価し、浮動小数点数に変換する。
詳細説明:
`eval()` 関数を使用しますが、`SAFE_GLOBALS` を指定することで、
実行可能な関数や変数を制限し、セキュリティを確保します。
例えば、"1/3" や "sqrt(2)" のような算術式を評価できます。
:param expr: str: 評価する数値式。
:returns: float: 評価された数値。
:raises ValueError: 無効な数値式が指定された場合。
"""
try:
return float(eval(expr, SAFE_GLOBALS, {}))
except Exception as e:
raise ValueError(f"Invalid numeric expression '{expr}': {e}")
[ドキュメント]
def parse_conversion_matrix(spec: str) -> Tuple[np.ndarray, np.ndarray]:
"""
変換行列と並進ベクトルを表す文字列を解析する。
詳細説明:
`'(a,b,c)(d,e,f)(g,h,i)'` または `'(a,b,c,tx)(d,e,f,ty)(g,h,i,tz)'`
形式の文字列を解析し、`numpy.ndarray` 形式の変換行列 `T` と
並進ベクトル `t` を返します。各要素は算術式として解釈されます。
:param spec: str: 変換行列と並進ベクトルを定義する文字列。
:returns: Tuple[np.ndarray, np.ndarray]: 変換行列 `T` と並進ベクトル `t` のタプル。
:raises ValueError: 不適切な形式の文字列が指定された場合。
"""
s = spec.replace(" ", "")
groups = re.findall(r'\(([^()]*)\)', s)
if len(groups) != 3:
raise ValueError("Matrix spec must have exactly 3 parenthesis groups.")
T = np.zeros((3,3), dtype=float)
t = np.zeros(3, dtype=float)
for i, g in enumerate(groups):
parts = g.split(',')
if len(parts) not in (3,4):
raise ValueError("Each row must have 3 or 4 entries.")
T[i,0] = parse_number(parts[0])
T[i,1] = parse_number(parts[1])
T[i,2] = parse_number(parts[2])
if len(parts) == 4:
t[i] = parse_number(parts[3])
return T, t
# ------------------------- Reporting -------------------------
[ドキュメント]
def site_effective_occupancy_and_mass_amu(site) -> Tuple[float, float]:
"""
サイトの有効占有率と合計質量を計算する。
詳細説明:
pymatgenのSiteオブジェクトから、原子種の占有率の合計（有効原子数）と
原子質量単位 (amu) での合計質量を計算します。
部分占有サイトも考慮されます。
:param site: pymatgen.core.structure.Site: 計算対象のサイトオブジェクト。
:returns: Tuple[float, float]: 有効占有率と合計質量 (amu) のタプル。
"""
occ = 0.0
mass = 0.0
for sp, amt in site.species.items():
occ += float(amt)
amu = getattr(sp, "atomic_mass", 0.0)
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]:
"""
構造全体の有効原子数、合計質量、サイト数を計算する。
詳細説明:
構造内の各サイトについて `site_effective_occupancy_and_mass_amu` を呼び出し、
全体の有効原子数 (occupancyの合計)、合計質量 (amu)、
およびサイトの総数を集計します。
:param struct: pymatgen.core.structure.Structure: 計算対象の構造オブジェクト。
:returns: Tuple[float, float, int]: 有効原子数、合計質量 (amu)、サイト数のタプル。
"""
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(struct: Structure, label: str) -> Tuple[float, float]:
"""
結晶構造の基本的な情報と密度を報告する。
詳細説明:
格子定数、格子角、セル体積、サイト数、有効原子数、総質量、
原子密度、質量密度を標準出力に表示します。
質量密度はamuからグラム、体積はÅ^3からcm^3に変換して計算されます。
:param struct: pymatgen.core.structure.Structure: 報告対象の構造オブジェクト。
:param label: str: 報告のヘッダーに使用されるラベル（例: "Original" または "Converted"）。
:returns: Tuple[float, float]: 原子密度 (atoms/Å^3) と質量密度 (g/cm^3) のタプル。
"""
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)
total_mass_g = total_mass_amu * 1.66053906660e-24 # g
volume_cm3 = volume * 1.0e-24 # cm^3
rho_atom = eff_atoms / volume if volume > 0 else float("nan")
rho_mass = total_mass_g / volume_cm3 if volume_cm3 > 0 else float("nan")
print(f"\n{label} Structure:")
print(f" Lattice constants : a={a:.6f} b={b:.6f} c={c:.6f}")
print(f" Lattice angles : alpha={alpha:.2f} beta={beta:.2f} gamma={gamma:.2f}")
print(f" Volume per cell : {volume:.6f} Å^3")
print(f" Sites (positions) : {n_sites}")
print(f" Effective atoms (Σ occ) : {eff_atoms:.6f}")
print(f" Total mass (cell) : {total_mass_amu:.6f} amu ({total_mass_g:.6e} g)")
print(f" Atomic density : {rho_atom:.6f} atoms/Å^3")
print(f" Mass density : {rho_mass:.6f} g/cm^3")
return rho_atom, rho_mass
[ドキュメント]
def dump_sites(struct: Structure, label: str):
"""
構造内の各サイトの情報を標準出力に表示する。
詳細説明:
各サイトのインデックス、原子種とその占有率、分数座標を整形して表示します。
:param struct: pymatgen.core.structure.Structure: サイト情報表示対象の構造オブジェクト。
:param label: str: 報告のヘッダーに使用されるラベル。
"""
print(f"\n--- {label} Structure Sites ---")
for i, site in enumerate(struct.sites):
sp_occ = ", ".join([f"{sp.symbol}:{amt:.3f}" for sp, amt in site.species.items()])
x, y, z = site.frac_coords
print(f"[{i:3d}] {sp_occ:20s} frac=({x:8.4f}, {y:8.4f}, {z:8.4f})")
# ------------------------- Symmetry helpers -------------------------
[ドキュメント]
def detect_setting_and_centering(struct: Structure, sym_tol: float) -> Tuple[str, str, str]:
"""
構造の空間群、R-格子設定、格子心タイプを検出する。
詳細説明:
`pymatgen.symmetry.analyzer.SpacegroupAnalyzer` を使用して空間群を特定し、
特にR-格子の場合の六方晶または菱面体晶設定、
および一般的な格子心タイプ (P, A, B, C, I, F, R) を識別します。
:param struct: pymatgen.core.structure.Structure: 分析対象の構造オブジェクト。
:param sym_tol: float: 対称性解析に使用する許容誤差。
:returns: Tuple[str, str, str]:
空間群シンボル、R-格子の設定（'hexagonal'または'rhombohedral'）、
格子心タイプ（'P', 'A', 'B', 'C', 'I', 'F', 'R'のいずれか）のタプル。
"""
sga = SpacegroupAnalyzer(struct, symprec=sym_tol)
spg = sga.get_space_group_symbol()
first = spg[0] if spg else '-'
setting = '-'
if first == 'R':
gamma = struct.lattice.gamma
setting = "hexagonal" if np.isclose(gamma, 120.0, atol=0.1) else "rhombohedral"
centering = first if first in {'P','A','B','C','I','F','R'} else '-'
return spg, setting, centering
[ドキュメント]
def hex_to_rhombo_T() -> np.ndarray:
"""
六方晶設定から菱面体晶設定への変換行列を返す。
詳細説明:
六方晶 (Hexagonal) 記述のR-格子を菱面体晶 (Rhombohedral) 記述の
原始セルに変換するための3x3変換行列を返します。
:returns: np.ndarray: 六方晶から菱面体晶への変換行列。
"""
return np.array([[ 2/3, 1/3, 1/3],
[-1/3, 1/3, 1/3],
[-1/3, -2/3, 1/3]], dtype=float)
[ドキュメント]
def rhombo_to_hex_T() -> np.ndarray:
"""
菱面体晶設定から六方晶設定への変換行列を返す。
詳細説明:
菱面体晶 (Rhombohedral) 記述のR-格子を六方晶 (Hexagonal) 記述の
慣用セルに変換するための3x3変換行列を返します。
:returns: np.ndarray: 菱面体晶から六方晶への変換行列。
"""
return np.array([[ 1, -1, 0],
[ 0, 1, -1],
[ 1, 1, 1]], dtype=float)
[ドキュメント]
def centering_to_prim_T(centering: str) -> Optional[np.ndarray]:
"""
格子心タイプに応じた原始セルへの変換行列を返す。
詳細説明:
面心 (F), 体心 (I), A/B/C面心 (A, B, C) の格子から
原始セルへ変換するための3x3変換行列を返します。
原始格子 (P) の場合は `None` を返します。
:param centering: str: 格子心タイプ（'F', 'I', 'A', 'B', 'C'など）。
:returns: Optional[np.ndarray]: 原始セルへの変換行列、または変換が不要な場合は `None`。
"""
c = centering.upper()
if c == 'F':
return np.array([[0.0, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.5, 0.5, 0.0]], dtype=float)
if c == 'I':
return np.array([[-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5]], dtype=float)
if c == 'A':
return np.array([[1.0, 0.0, 0.0],
[0.0, 0.5, 0.5],
[0.0,-0.5, 0.5]], dtype=float)
if c == 'B':
return np.array([[0.5, 0.0,-0.5],
[0.0, 1.0, 0.0],
[0.5, 0.0, 0.5]], dtype=float)
if c == 'C':
return np.array([[0.5,-0.5, 0.0],
[0.5, 0.5, 0.0],
[0.0, 0.0, 1.0]], dtype=float)
return None
# ------------------------- Basis change core -------------------------
[ドキュメント]
def change_basis_preserving_geometry(structure: Structure,
T: np.ndarray,
t: Optional[np.ndarray] = None,
xyz_tol: float = 1e-6) -> Structure:
"""
変換行列と並進ベクトルを用いて結晶構造の基底を変換する。
詳細説明:
新しい格子ベクトル `V'` = `T @ V` および、
新しい分数座標 `f'` = `f @ inv(T) + t` の変換を行います。
変換後のサイト座標は[0, 1)の範囲にラップされます。
変換後に同一とみなされる位置のサイトは結合され、重複が排除されます。
これにより、変換行列の行列式が1ではない場合でも、物理的な構造は維持されます。
:param structure: pymatgen.core.structure.Structure: 基底変換の対象となる元の構造。
:param T: np.ndarray: 新しい基底を定義する3x3の変換行列。
:param t: Optional[np.ndarray]: サイト座標に適用する並進ベクトル（デフォルトはなし）。
:param xyz_tol: float: 変換後に重複するサイトをマージするための座標の許容誤差。
:returns: pymatgen.core.structure.Structure: 基底変換とサイトのマージが適用された新しい構造。
:raises ValueError: 変換行列が特異な（行列式がほぼ0の）場合。
"""
print()
print("Transform lattice:")
detT = float(np.linalg.det(T))
if abs(detT) < 1e-12:
raise ValueError("Transformation matrix is singular (det ~ 0).")
V = np.array(structure.lattice.matrix, dtype=float)
Vp = T @ V
Tinv = np.linalg.inv(T)
frac = np.array([site.frac_coords for site in structure.sites], dtype=float)
frac_p = frac @ Tinv
if t is not None:
frac_p = frac_p + t.reshape(1,3)
frac_p = np.mod(frac_p, 1.0)
# 同一位置の原子のうち１つだけ残してkept_species/coordsに残す
round_ndp = max(0, int(-np.log10(max(xyz_tol, 1e-12)))) + 1
kept_coords: List[np.ndarray] = []
kept_species = []
seen = set()
for fc, site in zip(frac_p, structure.sites):
key = (round(fc[0]+xyz_tol, round_ndp), round(fc[1]+xyz_tol, round_ndp), round(fc[2]+xyz_tol, round_ndp),
site.species)
if key in seen: continue
seen.add(key)
kept_species.append(site.species)
kept_coords.append(fc)
new_lattice = Lattice(Vp)
new_structure = Structure(new_lattice, kept_species, kept_coords, coords_are_cartesian=False)
new_structure.sort()
removed = len(structure.sites) - len(new_structure.sites)
if removed > 0:
print(f"\n[Info] Merged {removed} duplicate sites due to det(T)={detT:.6f}.")
return new_structure
# ------------------------- Main -------------------------
[ドキュメント]
def main():
"""
メイン処理を実行し、結晶構造の変換と報告を行う。
詳細説明:
1. コマンドライン引数を解析し、入力CIFファイルを読み込む。
2. 元の構造の空間群、R-格子設定、格子心タイプを検出する。
3. 元の構造の情報を報告し、サイト情報を表示する。
4. 指定された変換タイプ（`prim`, `rhomb`, `hex`, `orth`, `MATRIX`）に応じて、
適切な変換行列を生成または解析する。
5. 変換行列と並進ベクトルを用いて構造の基底を変換する。
6. 変換後の構造情報を報告し、サイト情報を表示する。
7. 変換前後の原子密度と質量密度の整合性をチェックし、結果を標準出力に表示する。
8. 変換された構造を新しいCIFファイルとして保存する。
9. 部分占有に関する注意喚起メッセージを表示する。
"""
args = initialize()
infile = args.input_file
conv_spec = args.conversion
direction = args.direction
sym_tol = args.sym_tol
eps = args.eps
xyz_tol = args.xyz_tol
if not os.path.exists(infile):
print(f"Error: Input file '{infile}' not found.")
terminate()
print("Lattice conversion using pymatgen")
print(f"Input file : {infile}")
print(f"Conversion : {conv_spec}")
print(f"Direction : {direction}")
print(f"Symmetry tolerance : {sym_tol}")
print(f"(x,y,z) tolerance : {xyz_tol}")
print("-" * 30)
# Read CIF with occupancy_tolerance=0
try:
parser = CifParser(infile)
# parser = CifParser(infile, occupancy_tolerance = 1.0e-4) # 本来はここでpartial occupancyを許容すべきか否かを考慮
s_orig = parser.parse_structures(primitive=False)[0]
except Exception as e:
print(f"Error reading {infile}: {e}")
print(f" The CIF may include partial occupancy site.\n")
terminate()
spg, setting, centering = detect_setting_and_centering(s_orig, sym_tol)
print(f"Space Group: {spg}")
if setting != '-':
print(f"Detected setting : {setting} (R-lattice)")
print(f"Detected centering : {centering}")
rho_atom_orig, rho_mass_orig = report_structure(s_orig, "Original")
dump_sites(s_orig, "Original")
T = None
tvec = np.zeros(3, dtype=float)
conv_key = conv_spec.strip().lower()
if conv_key == 'prim':
print("\n[Conversion] primitive cell via SpacegroupAnalyzer")
try:
print("Transform lattice:")
s_conv = SpacegroupAnalyzer(s_orig, symprec=sym_tol).get_primitive_standard_structure()
except Exception:
# pymatgen 2022.9.21 以降では get_primitive_standard_structure() は pymatgen の内部変換に依存し、
# get_primitive_structure() は独自の実装を使う。
# 通常は get_primitive_standard_structure() で問題ないが、フォールバックとして get_primitive_structure() を用意。
s_conv = s_orig.get_primitive_structure()
rho_atom_new, rho_mass_new = report_structure(s_conv, "Converted")
dump_sites(s_conv, "Converted")
elif conv_key in ('rhomb', 'hex', 'orth'):
print("\n[Conversion] auto rule:", conv_key)
if conv_key == 'rhomb' and centering == 'R' and setting == 'hexagonal':
T = hex_to_rhombo_T()
elif conv_key == 'hex' and centering == 'R' and setting == 'rhombohedral':
T = rhombo_to_hex_T()
elif conv_key == 'orth':
T = centering_to_prim_T(centering)
else:
T = np.identity(3) # 適用できる変換がない場合は恒等変換
else:
print("\n[Conversion] explicit matrix parsing")
T, tvec = parse_conversion_matrix(conv_spec)
if direction == "ConvertedToOriginal":
T = np.linalg.inv(T)
tvec = -tvec @ T
if conv_key != 'prim' and T is None:
print()
print(f"Error in main(): Inconsistent lattice for concersion to [{conv_key}]")
print()
terminate()
# prim 以外の変換、またはカスタム行列が指定された場合に change_basis_preserving_geometry を実行
if conv_key != 'prim' and T is not None:
print_matrix("Transformation matrix [T]:", T)
print(f"det(T) = {np.linalg.det(T):.8f}")
s_conv = change_basis_preserving_geometry(s_orig, T, t=tvec, xyz_tol=xyz_tol)
rho_atom_new, rho_mass_new = report_structure(s_conv, "Converted")
dump_sites(s_conv, "Converted")
print()
print("Check consistency:")
# 密度比較は相対誤差でチェック
if abs((rho_atom_new - rho_atom_orig) / rho_atom_orig) < eps:
print(f" Atomic densities are identical ({rho_atom_orig:.6f} vs. {rho_atom_new:.6f} atoms/Å^3) within eps={eps}")
else:
print("#" * 80)
print(f" Error!!: Atomic density changed from {rho_atom_orig:.6f} to {rho_atom_new:.6f} atoms/Å^3")
print("#" * 80)
print()
if abs((rho_mass_new - rho_mass_orig) / rho_mass_orig) < eps:
print(f" Mass densities are identical ({rho_mass_orig:.6f} vs. {rho_mass_new:.6f} g/cm^3) within eps={eps}")
else:
print("#" * 80)
print(f" Error!!: Mass density changed from {rho_mass_orig:.6f} to {rho_mass_new:.6f} g/cm^3")
print("#" * 80)
print()
outfile = os.path.splitext(infile)[0] + "_converted.cif"
s_conv.to(filename=outfile, fmt="cif")
print(f"\nConverted CIF file saved to: {outfile}")
print()
print("#"*80)
print("# Please check OCCUPANCY: pymatgen may fail to read partial occupancies ocrrectly.")
print("#"*80)
if __name__ == "__main__":
main()
input("\nPress ENTER to terminate>>")