cif2rdf.py 技術ドキュメント

cif2rdf.py は、CIF (Crystallographic Information File) 形式の結晶構造データから動径分布関数 (RDF) と配位数 (CN) を計算し、その結果をCSVファイルとして保存し、グラフとして可視化するPythonプログラムです。

プログラムの動作

このプログラムは、入力されたCIFファイルから結晶構造情報を読み取り、指定されたパラメータに基づいて動径分布関数 (RDF) および累積配位数 (RCN) を計算します。計算されたデータはCSVファイル形式で出力され、同時にMatplotlibライブラリを使用してグラフとして可視化されます。

主な機能は以下の通りです。

  • CIFファイルの解析: 結晶構造データを含むCIFファイルを読み込み、格子定数、原子サイト、原子種などの情報を抽出します。

  • 動径分布関数 (RDF) の計算: 結晶構造内の指定された原子またはサイトを中心として、他の原子までの距離分布を計算します。周期境界条件が考慮されます。

  • 累積配位数 (RCN) の計算: RDFを積分し、指定された距離内にある原子の総数を計算します。

  • 計算モードの選択:

    • atom モード: 各原子種 (例: O, Ga) ごとのRDFとRCNを計算します。

    • site モード: 各非対称原子サイト (例: O1, Ga1) ごとのRDFとRCNを計算します。

    • cn モード: 各原子種ごとのRDFとRCNに加え、配位数ごとのCNRDFとCNRCNを計算します。

  • ガウシアン畳み込み: オプションとして、計算されたRDFに対してガウシアン関数による畳み込み（スムージング）を適用できます。

  • データ出力: 計算結果をCSVファイルとして保存します。

  • グラフ表示: 計算されたRDFおよびRCNのグラフをリアルタイムで表示します。

このプログラムは、結晶材料における原子の局所的な環境や配位構造を解析し、構造的な特徴を理解するためのツールとして機能します。

原理

動径分布関数 (RDF)

プログラムでは、まず結晶構造内の各非対称原子サイト (atom0) を中心として、他の全ての原子 (atom1) までの距離を計算します。この距離計算では、結晶の周期性を考慮するために、隣接する単位胞内の原子も探索します。

動径分布関数 \(RDF_{i \to j}(R)\) は、中心原子 \(i\) から距離 \(R\) の球面殻内に存在する原子種 \(j\) の数密度に比例する量を示します。本プログラムでは、指定されたRの範囲 \(R_{min}\) から \(R_{max}\) までを微小なステップ \(R_{step}\) で区切り、各距離ビンに属する原子の数をカウントすることでRDFを計算します。

具体的には、中心原子 atom0 から検出された原子 atom1 までの距離 dis が \(R_{min} \leq dis \leq R_{max}\) の範囲にある場合、その距離 dis に対応するビン idx = int(dis / Rstep) に atom1 の占有率 occ1 を加算します。

最終的に、各ビンの値は \(R_{step}\) で割って正規化されます。 $\(RDF_{type}(R) = \frac{\sum_{i \in type} N_{i}(R)}{R_{step}}\)\( ここで、\)N_i(R)\( は中心原子 \)i\( から距離 \)R$ にある原子のカウントです。この正規化により、単位距離あたりの原子数を表現します。

累積配位数 (RCN)

累積配位数 (RCN) は、動径分布関数を積分したもので、ある距離 \(R\) 以内に存在する原子の総数を示します。プログラムでは、RDFの離散的な合計としてRCNを計算します。 $\(RCN(R) = \sum_{r' \le R} RDF(r') \cdot R_{step}\)\( ここで、\)RDF(r')\( は距離 \)r'$ における動径分布関数の値です。

配位数 (CN) 別RDF/RCN (mode='cn')

mode='cn' の場合、特定の配位数を持つ原子に限定したRDFおよびRCNを計算します。これは、各中心原子についてRCNを計算し、そのRCNが設定された閾値（例: \(CN \ge 1\), \(CN \ge 2\), ..., \(CN \ge nCN\)）を超える距離において、その原子がその配位数以上に寄与するものとしてカウントすることで実現されます。

ガウシアン畳み込み

計算されたRDFは、ノイズを低減し、より滑らかな曲線を得るために、ガウシアン関数による畳み込み処理が適用されることがあります (wConv が0でない場合)。ガウシアン関数 \(G(x, \mu, \sigma)\) は以下の形で与えられます。 $\(G(x, \mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\)$ 畳み込み処理は tklib.tksci.tkconvolution.convolute_by_func 関数によって行われます。ここで wConv はガウシアン関数の半値幅に相当するパラメータとして使用されます。

必要な非標準ライブラリとインストール方法

このプログラムは以下の非標準ライブラリに依存しています。

  • numpy: 高度な数値計算（配列操作など）のための基盤ライブラリ。

  • scipy: 科学技術計算ライブラリ。特に補間機能 scipy.interpolate.interp1d が使用されています。

  • matplotlib: グラフ描画のためのライブラリ。

  • tklib: このプログラムの作者によって提供されるカスタムライブラリ。ファイル操作 (tkFile)、ユーティリティ関数 (tkutils)、科学計算用モジュール (tkmatrix, tkconvolution)、および結晶学関連のクラス (tkCIF, tkCrystal, tkAtomType) を含みます。

numpy, scipy, matplotlib は Python Package Index (PyPI) に登録されており、pip コマンドを使ってインストールできます。

pip install numpy scipy matplotlib

tklib はPyPIに登録されていないカスタムライブラリであるため、pip コマンドによる直接的なインストールはできません。このライブラリは、cif2rdf.py プログラムと同じ環境に別途提供されるか、適切にPythonのパスが設定された場所に配置されている必要があります。

必要な入力ファイル

プログラムは以下の入力ファイルを必要とします。

  • CIFファイル (.cif):

    • 形式: Crystallographic Information File (CIF) 形式。国際結晶学連合 (IUCr) によって標準化されたテキストファイル形式で、結晶構造に関する詳細情報（格子定数、空間群、原子座標、原子種、占有率など）を含みます。

    • ファイル名: コマンドライン引数で指定します。指定がない場合、プログラム内のデフォルト値 (IGZO.cif) が使用されます。

    • 内容: tklib.tkcrystal.tkcif.tkCIF クラスによって解析され、結晶構造のオブジェクト (tklib.tkcrystal.tkcrystal.tkCrystal) が構築されます。

生成される出力ファイル

プログラムは、計算結果を以下の形式で出力します。

  1. RDF/RCNデータCSVファイル: [入力CIFファイル名から拡張子を除いた部分]-RDF.csv

    • 例: IGZO.cif を入力した場合、IGZO-RDF.csv

    • 内容:

      • 1列目: R(A) (距離、オングストローム単位)。

      • 続く列: 計算モードによって異なります。

        • mode='atom' または mode='cn' の場合: 各原子種 (例: O, Ga, In, Zn) ごとのRDF値とRCN値。列名は RDF(AtomTypeName:index) および RCN(AtomTypeName:index)。

        • mode='site' の場合: 各非対称原子サイト (例: O1, Ga1) ごとのRDF値とRCN値。列名は RDF(AtomSiteName:index) および RCN(AtomSiteName:index)。

  2. 配位数別CNRDF/CNRCNデータCSVファイル (モードがcnの場合のみ):

    • [入力CIFファイル名から拡張子を除いた部分]-CNRDF-[AtomTypeName].csv

    • 例: IGZO.cif を入力し、原子種が O の場合、IGZO-CNRDF-O.csv

    • 内容:

      • 1列目: R(A) (距離、オングストローム単位)。

      • 続く列: 各原子種について、設定された配位数 (nCN) までの各閾値 (icn = 1, ..., nCN) におけるCNRDF値とCNRCN値。列名は CNRDF(CN>=k)(AtomTypeName:index) および CNRCN(CN>=k)(AtomTypeName:index)。

  3. グラフウィンドウ:

    • 計算されたRDFとRCNのグラフがMatplotlibによって表示されます。

    • mode='atom' または mode='site' の場合、RDFとRCNがそれぞれ1つのプロットにまとめられます。

    • mode='cn' の場合、全体のRDF/RCNと、各原子種ごとの配位数別CNRDF/CNRCNが個別のサブプロットに表示されます。このウィンドウを閉じるとプログラムは終了します。

コマンドラインでの使用例 (Usage)

プログラムはコマンドライン引数を受け付けて動作を制御します。

Usage:
  (i) python cif2rdf.py mode ciffile Rmin Rmax nR nCN wConv
      mode: [atom|site|cn]
      ex: python cif2rdf.py cn IGZO.cif 0.1 5.0 1001 6 0.0

  • mode: 計算モードを指定します。以下のいずれかを選択します。

    • atom: 原子種ごとのRDF/RCNを計算します。

    • site: 非対称原子サイトごとのRDF/RCNを計算します。

    • cn: 原子種ごとのRDF/RCNに加え、配位数別CNRDF/CNRCNを計算します。

  • ciffile: 解析するCIFファイルのパス。

  • Rmin: RDF/RCN計算の開始距離 (Å)。

  • Rmax: RDF/RCN計算の終了距離 (Å)。

  • nR: \(R_{min}\) から \(R_{max}\) までの距離を分割する点数。Rstep は Rmax / (nR - 1) で計算されます。

  • nCN: mode='cn' の場合に計算する最大配位数。1 から nCN までの配位数について結果が出力されます。他のモードではこの値は無視されますが、コマンドライン引数としては任意の整数（例: 0 や 1）を指定する必要があります。

  • wConv: ガウシアン畳み込みの半値幅 (Å)。0.0 を指定すると畳み込みは行われません。

コマンドラインでの具体的な使用例

ここでは、IGZO.cif というCIFファイルが存在すると仮定し、mode='cn' で、Rの範囲を0.1 Åから5.0 Å、点数を1001点、最大配位数を6、畳み込みなしで実行する例を示します。

python cif2rdf.py cn IGZO.cif 0.1 5.0 1001 6 0.0

実行結果の説明

このコマンドを実行すると、以下の処理が行われます。

  1. プログラムは IGZO.cif ファイルを読み込み、結晶構造情報を解析します。

  2. 各原子種（例: O, Ga, In, Znなど、IGZO.cif 内の原子種による）について、RminからRmaxまでのRDFとRCNを計算します。

  3. さらに、各原子種について、配位数1から6までのCNRDFとCNRCNが計算されます。

  4. 計算結果は以下のCSVファイルに保存されます。

    • IGZO-RDF.csv: 全体のRDFおよびRCNデータ。

    • IGZO-CNRDF-O.csv: O原子に関する配位数別CNRDF/CNRCNデータ。

    • IGZO-CNRDF-Ga.csv: Ga原子に関する配位数別CNRDF/CNRCNデータ。

    • (他の原子種についても同様にファイルが生成されます)

  5. Matplotlibによって、計算されたRDF、RCN、配位数別CNRDF、CNRCNをプロットしたグラフウィンドウが表示されます。このウィンドウを閉じるとプログラムは終了します。

これにより、IGZO結晶構造における各原子種の局所的な構造情報が定量的に評価され、視覚的に確認できるようになります。