プログラム convert_cell.py 技術ドキュメント
プログラムの動作
convert_cell.py は、Pythonの pymatgen ライブラリを使用して、結晶構造記述ファイル（CIF形式）で定義された単位胞を別の格子設定に変換するためのプログラムです。このプログラムは、部分占有率を持つ原子サイトも適切に処理し、変換前後で物理的な密度が維持されていることを検証します。
主な機能は以下の通りです。
CIFファイルの読み込み: 入力として与えられたCIFファイルを解析し、結晶構造データ（格子定数、原子の種類、位置、占有率など）を読み込みます。
空間群解析と単位胞の検出:
pymatgen.symmetry.analyzer.SpacegroupAnalyzerを用いて、元の構造の空間群、R-格子の設定（六方晶または菱面体晶）、および中心化の種類（例: P, A, B, C, I, F）を自動的に検出します。格子変換の実行: 以下のいずれかの方法で単位胞を変換します。
prim: 空間群解析に基づいてプリミティブセル（原始単位胞）に変換します。rhomb: R-格子の六方晶（hexagonal）設定の慣用セルを、菱面体晶（rhombohedral）プリミティブセルに変換します。hex: R-格子の菱面体晶設定のプリミティブセルを、六方晶の慣用セルに変換します。orth: 体心（I）、面心（F）、底心（A, B, C）などの中心化された格子を、対応するプリミティブセルに変換します。MATRIX: ユーザーが任意の3x3または3x4の変換行列と並進ベクトルを指定して変換します。行列の各要素には、1/3やsqrt(2)/2のような算術式を使用できます。
構造情報の報告: 変換前後の構造について、以下の詳細情報を出力します。
格子定数 (\(a, b, c, \alpha, \beta, \gamma\))
単位胞の体積 (Å\(^3\))
サイト数（位置の数）
有効原子数（原子占有率の合計）
単位胞の総質量 (amu, g)
原子密度 (atoms/Å\(^3\))
質量密度 (g/cm\(^3\))
密度の整合性チェック: 変換前後で計算された原子密度と質量密度が、指定された相対許容誤差
eps(デフォルト 1e-4) の範囲内で一致することを確認します。これにより、変換が物理的に妥当であるかを検証します。一致しない場合は警告メッセージを出力します。変換後のCIFファイルの出力: 変換された構造を新しいCIFファイルとして保存します。
このプログラムは、結晶構造シミュレーションや解析において、特定の格子設定への変換が必要な場面、あるいは異なる格子設定の構造データを比較する際に役立ちます。
原理
convert_cell.py は、結晶学における単位胞の基底変換と、それに伴う原子座標の変換を数学的に実行します。
1. 基底ベクトルと変換行列
結晶構造の単位胞は、3つの基底ベクトル \(\mathbf{a}, \mathbf{b}, \mathbf{c}\) によって定義されます。これらをPythonの numpy 配列として扱う際、pymatgen.core.lattice.Lattice オブジェクトの matrix プロパティは、これらのベクトルを行ベクトルとして持つ3x3行列 \(V\) を返します。
ここで、\(\mathbf{a}, \mathbf{b}, \mathbf{c}\) はそれぞれ \( (a_x, a_y, a_z) \), \( (b_x, b_y, b_z) \), \( (c_x, c_y, c_z) \) といったデカルト座標成分を持つ行ベクトルです。
新しい単位胞の基底ベクトル \(\mathbf{a}', \mathbf{b}', \mathbf{c}'\) は、変換行列 \(T\) (3x3行列) を用いて元の基底ベクトルの線形結合として表現されます。これにより、新しい基底ベクトルを行ベクトルとして持つ行列 \(V'\) は以下のように計算されます。
具体的には、各新しい基底ベクトルは次のように定義されます。 $\( \mathbf{a}' = T_{11}\mathbf{a} + T_{12}\mathbf{b} + T_{13}\mathbf{c} \\ \mathbf{b}' = T_{21}\mathbf{a} + T_{22}\mathbf{b} + T_{23}\mathbf{c} \\ \mathbf{c}' = T_{31}\mathbf{a} + T_{32}\mathbf{b} + T_{33}\mathbf{c} \)$
2. 原子分数座標の変換
新しい単位胞 \(V'\) 内における原子の分数座標 \(\mathbf{f}' = (x', y', z')\) は、元の単位胞 \(V\) 内における分数座標 \(\mathbf{f} = (x, y, z)\)、変換行列 \(T\) の逆行列 \(T^{-1}\)、および並進ベクトル \(\mathbf{t} = (t_x, t_y, t_z)\) を用いて計算されます。
ここで \(\mathbf{f}\) は行ベクトルとして扱われます。変換後、分数座標は周期境界条件に従って \(0 \le x, y, z < 1\) の範囲に正規化されます（np.mod 関数を使用）。
変換行列 \(T\) の行列式 \( \det(T) \) は、単位胞の体積比に関係します。新しい単位胞の体積 \(V_{new}\) と元の単位胞の体積 \(V_{old}\) の間には以下の関係があります。 $\( V_{new} = |\det(T)| V_{old} \)$
また、変換によって複数の原子サイトが新しい単位胞内の同一の分数座標にマッピングされる場合（特に \(|\det(T)| < 1\) のとき）、それらの重複サイトは xyz_tol で指定された許容誤差内で一つに統合されます。
3. 密度計算
プログラムは、変換前後の構造について以下の密度関連量を計算し、物理的な整合性を確認します。
有効原子数: 各サイトにおける原子の占有率（occupancy）の合計として計算されます。部分占有率を持つ原子は、その占有率分だけ原子数に寄与します。 $\( N_{eff} = \sum_{i \in \text{sites}} \sum_{sp \in \text{species}_i} \text{occupancy}(sp) \)$
単位胞の総質量: 各原子サイトの元素の原子質量（amu）と占有率を乗じて合計することで計算されます。その後、グラム単位に変換されます。 $\( M_{total, amu} = \sum_{i \in \text{sites}} \sum_{sp \in \text{species}_i} \text{occupancy}(sp) \cdot \text{atomic\_mass}(sp) \)\( ここで、\)1 \text{ amu} = 1.66053906660 \times 10^{-24} \text{ g}$ を用いてグラム単位に変換されます。
体積:
pymatgenが計算するÅ\(^3\)単位の体積を、質量密度計算のために \(\text{cm}^3\) に変換します。 \(1 \text{ Å}^3 = 1.0 \times 10^{-24} \text{ cm}^3\)原子密度: $\( \rho_{atom} = \frac{N_{eff}}{\text{体積}(\text{Å}^3)} \)$
質量密度: $\( \rho_{mass} = \frac{M_{total}(\text{g})}{\text{体積}(\text{cm}^3)} \)$
4. 組み込み変換行列
特定の変換オプション (rhomb, hex, orth) に対して、プログラムは以下の固定変換行列を使用します。
Hexagonal to Rhombohedral (R-格子の六方晶設定から菱面体晶プリミティブセルへ): $\( T_{Hex \to Rhomb} = \begin{pmatrix} 2/3 & 1/3 & 1/3 \\ -1/3 & 1/3 & 1/3 \\ -1/3 & -2/3 & 1/3 \end{pmatrix} \)$
Rhombohedral to Hexagonal (R-格子の菱面体晶設定から六方晶慣用セルへ): $\( T_{Rhomb \to Hex} = \begin{pmatrix} 1 & -1 & 0 \\ 0 & 1 & -1 \\ 1 & 1 & 1 \end{pmatrix} \)$
Centered to Primitive (面心、体心、底心格子からプリミティブセルへ):
F-centering (面心): $\( T_F = \begin{pmatrix} 0 & 1/2 & 1/2 \\ 1/2 & 0 & 1/2 \\ 1/2 & 1/2 & 0 \end{pmatrix} \)$
I-centering (体心): $\( T_I = \begin{pmatrix} -1/2 & 1/2 & 1/2 \\ 1/2 & -1/2 & 1/2 \\ 1/2 & 1/2 & -1/2 \end{pmatrix} \)$
A-centering (A-底心): $\( T_A = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1/2 & 1/2 \\ 0 & -1/2 & 1/2 \end{pmatrix} \)$
B-centering (B-底心): $\( T_B = \begin{pmatrix} 1/2 & 0 & -1/2 \\ 0 & 1 & 0 \\ 1/2 & 0 & 1/2 \end{pmatrix} \)$
C-centering (C-底心): $\( T_C = \begin{pmatrix} 1/2 & -1/2 & 0 \\ 1/2 & 1/2 & 0 \\ 0 & 0 & 1 \end{pmatrix} \)$
これらの変換は、結晶学の標準的な手法に基づいており、pymatgen の SpacegroupAnalyzer も同様の原理でプリミティブセルを導出します。MATRIX オプションでは、ユーザーが任意の変換行列を定義できるため、より柔軟な変換が可能です。
必要な非標準ライブラリとインストール方法
convert_cell.py の実行には、以下のPython非標準ライブラリが必要です。
pymatgen: 結晶構造の解析、CIFファイルの読み書き、格子変換のコア機能を提供します。numpy: 数値計算、特に多次元配列操作や線形代数演算（行列積、逆行列計算など）に使用されます。
これらのライブラリは、Pythonのパッケージマネージャーである pip を使ってインストールできます。
pip install pymatgen numpy
必要な入力ファイル
ファイル名: コマンドライン引数で指定される任意のCIFファイルパス。
形式: CIF (Crystallographic Information File) 形式。
ファイルは、少なくとも以下の結晶構造情報を含んでいる必要があります。
格子定数 (a, b, c, alpha, beta, gamma)
原子の種類（元素記号）
原子の分数座標 (x, y, z)
原子の占有率（オプションですが、部分占有率を扱う場合は必須）
プログラムは
pymatgen.io.cif.CifParserを使用してファイルを読み込みます。このパーサーは、デフォルトでoccupancy_tolerance=0を設定しているため、厳密な占有率の値を持つサイトも適切に解釈します。
生成される出力ファイル
ファイル名: 入力CIFファイルの名前の末尾に
_converted.cifを追加した名前で保存されます。例: 入力ファイルが
structure.cifの場合、出力ファイルはstructure_converted.cifとなります。
内容: 変換後の結晶構造情報が記述されたCIFファイルです。
新しい格子定数と格子角度。
変換後の単位胞における原子サイトの分数座標。
原子の種類とそれぞれの占有率。
結晶の対称性情報（空間群など）は、変換方法によっては変更される可能性があります。
コマンドラインでの使用例 (Usage)
convert_cell.py は、以下の基本的な構文でコマンドラインから実行します。
python convert_cell.py <input_file> [-c {prim|rhomb|hex|orth|MATRIX}] [-d {OriginalToConverted|ConvertedToOriginal}] [-t SYMTOL] [-x XYZTOL] [--eps EPS]
引数の説明:
<input_file>:入力とするCIFファイルへのパス（必須）。
-c,--conversion:実行する格子変換の種類を指定します。
prim: 入力構造のプリミティブセルを生成します。SpacegroupAnalyzerを使用します。rhomb: R-格子の六方晶（hexagonal）設定を菱面体晶（rhombohedral）プリミティブセルに変換します。hex: R-格子の菱面体晶（rhombohedral）設定を六方晶（hexagonal）慣用セルに変換します。orth: 結晶の中心化（A, B, C, F, I）に基づいてプリミティブセルに変換します。MATRIX: ユーザー定義の変換行列を指定します。以下の2つの形式があります。'(a,b,c)(d,e,f)(g,h,i)': 3x3の変換行列。'(a,b,c,tx)(d,e,f,ty)(g,h,i,tz)': 3x3の変換行列と並進ベクトル \((t_x, t_y, t_z)\)。 各a,b,c,...は、1/2,sqrt(2)/2,piなどの算術式として記述できます。
-d,--direction:--conversion MATRIXオプションを使用する際に、変換行列の適用方向を指定します。OriginalToConverted(デフォルト): 指定されたMATRIXを直接、元の構造に適用します。ConvertedToOriginal: 指定されたMATRIXの逆行列を元の構造に適用します。
-t,--sym_tol:空間群解析を行う際の対称性許容誤差をÅ単位で指定します。デフォルト値は
0.01です。
-x,--xyz_tol:格子変換後に、同一位置（分数座標）に存在する重複サイトを統合する際の許容誤差を指定します。デフォルト値は
1.0e-4です。
--eps:変換前後で原子密度および質量密度を比較する際の相対許容誤差を指定します。デフォルト値は
1.0e-4です。
コマンドラインでの具体的な使用例
ここでは、いくつかの具体的な使用例とその実行結果について説明します。
例1: プリミティブセルへの変換
example.cif という名前の入力CIFファイルが存在すると仮定し、この構造をプリミティブセルに変換します。
python convert_cell.py example.cif -c prim
実行結果の説明:
プログラムは
example.cifを読み込み、pymatgenのSpacegroupAnalyzerを用いて、元の構造の空間群、R-格子設定（もしR-格子の場合）、および中心化の種類を検出して表示します。"Original Structure" として、元の単位胞の格子定数、体積、サイト数、有効原子数、総質量、原子密度、質量密度といった詳細な情報が報告されます。また、各原子サイトの分数座標も表示されます。
SpacegroupAnalyzerによってプリミティブセルへの変換が実行されます。"Converted Structure" として、変換後のプリミティブセルの同様の詳細情報が報告されます。元の構造に比べて体積が小さく、サイト数も減っていることが一般的です。
変換前後で原子密度と質量密度が
--eps(デフォルト 1e-4) 以内で一致しているかどうかが検証され、結果が報告されます。物理的な量が維持されていることを確認できます。最終的に、変換後の構造が
example_converted.cifという名前で新しいCIFファイルとして保存されます。最後に、部分占有率の読み取りに関する注意喚起メッセージが表示されます。
例2: ユーザー定義行列による変換 (3x3行列)
入力ファイル example.cif に対して、例えば \(c\) 軸の長さを半分にする変換行列（これにより体積が半分になる）を適用します。
python convert_cell.py example.cif -c '(1,0,0)(0,1,0)(0,0,1/2)'
実行結果の説明:
example.cifの元の構造情報が表示されます。指定された変換行列 \(T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0.5 \end{pmatrix}\) が解析され、これが元の構造に適用されます。
det(T)は0.5となり、新しい単位胞の体積は元の体積の半分になります。"Converted Structure" として、新しい格子定数、体積（元の約半分）、原子数、および密度が報告されます。
密度の一貫性チェックでは、体積が変化しても、単位体積あたりの原子数や質量は変わらないため、原子密度と質量密度は元の構造と一致することが確認されます。
変換後の構造は
example_converted.cifとして出力されます。
例3: ユーザー定義行列による変換と並進 (3x4行列)
入力ファイル example.cif に対して、単位胞の形状は変えずに原子の分数座標を全体的に \((0.1, 0.2, 0.3)\) だけ並進させる変換を適用します。
python convert_cell.py example.cif -c '(1,0,0,0.1)(0,1,0,0.2)(0,0,1,0.3)'
実行結果の説明:
example.cifの元の構造情報が表示されます。指定された変換行列 \(T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}\) と並進ベクトル \(\mathbf{t} = (0.1, 0.2, 0.3)\) が適用されます。
det(T)は1.0となり、体積は変化しません。"Converted Structure" では、格子定数や体積は変化せず、原子の分数座標が指定された並進ベクトル分だけシフトします。これに伴い、原子密度と質量密度も元の構造と完全に一致します。
変換後の構造は
example_converted.cifとして出力されます。
例4: R-格子の六方晶設定を菱面体晶プリミティブセルへ変換
r_hex.cif という名前のR-格子の六方晶設定を持つCIFファイルが存在すると仮定します。この構造を菱面体晶のプリミティブセルに変換します。
python convert_cell.py r_hex.cif -c rhomb
実行結果の説明:
r_hex.cifが読み込まれ、プログラムはR-格子の六方晶設定であることを検出します。元の構造情報が表示された後、Hexagonal to Rhombohedral 変換行列が自動的に適用されます。
"Converted Structure" として、菱面体晶プリミティブセルの格子定数、体積、原子数、および密度が報告されます。この変換により、元の六方晶セルよりも体積が小さく、より高い対称性を示すプリミティブセルが得られます。
密度の整合性が確認され、
r_hex_converted.cifが出力されます。