make_md_history.py 技術ドキュメント

プログラムの動作

make_md_history.py は、VASP (Vienna Ab-initio Simulation Package) の分子動力学 (MD) 計算の出力ファイルを解析し、主に以下の機能を提供します。

XCrySDen AXSF (アニメーションXSF) ファイルの生成: VASPの OUTCAR ファイルに含まれる各MDステップの結晶構造情報を抽出し、XCrySDenで可視化可能なアニメーション形式のAXSFファイル (md.xsf) を作成します。このファイルはMDシミュレーション中の原子の動きを視覚的に追跡するのに役立ちます。また、同時にMDステップごとの温度、エネルギーなどの情報が記載されたCSVファイル (inf.csv) も出力します。
平均二乗変位 (MSD) および関連履歴の計算: AXSFファイル (または必要に応じて OUTCAR から再生成) を使用して、MDシミュレーション中の原子の平均二乗変位 (MSD) を計算します。MSDは原子の拡散挙動を評価するための重要な指標です。MSDの計算に加えて、重心の変位、格子定数、体積、温度、エネルギーなどの時系列データをまとめた履歴CSVファイル (sample_name-history.csv) を生成します。初期構造 (POSCAR) と最終構造 (CONTCAR) はCIF形式でも保存されます。

このプログラムは、VASPのMDシミュレーション結果の解析と可視化を効率化し、材料の動的特性の理解を深めることを目的としています。

原理

1. 結晶構造の抽出とAXSF形式への変換

VASPの OUTCAR ファイルは、各MDステップにおける格子ベクトル、原子位置、原子にかかる力（または速度）などの情報を含んでいます。make_md_history.py は、これらの情報を逐次読み込み、tkCrystal オブジェクトとしてメモリ上に保持します。 AXSFファイルは、XCrySDenという結晶構造可視化ソフトウェアがサポートするアニメーション形式です。各MDステップの結晶構造を、以下の形式で記述します。

ANIMSTEPS [総ステップ数]
CRYSTAL
PRIMVEC [ステップ番号]
 [格子ベクトル a]
 [格子ベクトル b]
 [格子ベクトル c]
CONVVEC [ステップ番号]
 [格子ベクトル a]
 [格子ベクトル b]
 [格子ベクトル c]
PRIMCOORD [ステップ番号]
 [総原子数] 1
 [原子番号] [x座標] [y座標] [z座標] [Fx] [Fy] [Fz]
 ...

プログラムは OUTCAR から読み込んだ格子情報と原子座標（カルテシアン座標に変換）をこの形式に変換し、連続的にファイルに書き込むことでアニメーションを生成します。

2. 平均二乗変位 (MSD) の計算

MSDは、MDシミュレーションにおける原子の拡散特性を評価するための中心的な量です。時刻 $t$ における原子 $i$ の位置 $\mathbf{r}_i(t)$ と、基準時刻 $t_0$ における初期位置 $\mathbf{r}_i(t_0)$ を用いて、以下の式で定義されます。

\[ \text{MSD}(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} \| \mathbf{r}_i(t) - \mathbf{r}_i(t_0) \|^2 \right\rangle \]

ここで、$N$ は平均を取る原子の総数、$\| \cdot \|$ はユークリッド距離を示します。本プログラムの実装では、アンサンブル平均 ( $\left\langle \cdot \right\rangle$ ) は陽に計算されず、単に全MDステップにおける対象原子の変位の二乗和を $N$ で割る形を取ります。

a. 周期境界条件の取り扱いとオフセット補正

MDシミュレーションでは周期境界条件が適用されるため、原子は結晶の境界を越えて移動すると、反対側から再出現します。このため、単純に現在の座標と初期座標の差を取ると、実際の移動距離が大幅に過大評価される可能性があります。本プログラムでは、これを解決するために オフセット補正 を行います。各MDステップにおいて、原子が周期境界を越えたかどうかを検出し、その原子の「真の」絶対座標を追跡するために、仮想的な格子単位数分のオフセットを累積加算します。具体的には、ある原子が以前のステップから現在のステップに移動した際に、周期境界内での最短距離だけでなく、境界を越えた場合の変位（オフセット）を計算し、累積していきます。

b. 重心移動補正 (`drift_correction`)

シミュレーション全体で系が一定方向に移動（ドリフト）する場合があります。このドリフトは原子のランダムな拡散運動とは異なるため、MSD計算に影響を与えます。drift_correction オプションが 'cm' (Center of Mass) に設定されている場合、プログラムは各MDステップにおける系の重心を計算し、その重心の移動を各原子の変位から差し引くことで、ドリフトの影響を補正します。

系の重心 $\mathbf{R}_{\text{CM}}(t)$ は次式で計算されます。 $$ \mathbf{R}_{\text{CM}}(t) = \frac{\sum_i m_i \mathbf{r}_i(t)}{\sum_i m_i} $$ ここで、$m_i$ は原子 $i$ の質量です。補正された原子位置 $\mathbf{r}'_i(t)$ は、 $$ \mathbf{r}'_i(t) = \mathbf{r}_i(t) - (\mathbf{R}_{\text{CM}}(t) - \mathbf{R}_{\text{CM}}(t_0)) $$ となります。ただし、本プログラムの実装では、rglist[istep][i] が重心の移動量を示しているため、重心の絶対座標ではなく、pos[i] - offset[istep][isite][i] - rglist[istep][i] のように各原子のオフセット補正された位置から重心の位置そのものを減算しています。これにより、系の重心が原点に固定されたかのような相対座標でMSDが計算されます。

c. 基準位置 (`origin`)

MSDの計算における初期位置 $\mathbf{r}_i(t_0)$ は、origin オプションによって選択可能です。

'POSCAR' : VASPの初期構造ファイル POSCAR の原子位置を基準とします。
'OUTCAR' : VASP OUTCAR ファイルから読み込んだ最初のMDステップの原子位置を基準とします。
'average' : シミュレーションの最初の nskip ステップをスキップし、それ以降の原子位置を平均したものを基準とします。これは、初期緩和段階の非平衡状態を避け、より平衡に近い状態からの拡散を評価する際に有用です。

d. MSDの集計モード (`mode`)

MSDの集計方法を mode オプションで選択できます。

'atom' : 同じ原子種に属する全ての原子のMSDを平均します。
'site' : 個々の原子サイトごとにMSDを計算します。これにより、特定の原子がどの程度拡散したかを個別に追跡できます。

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

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

numpy: 数値計算ライブラリ。
```
pip install numpy
```

また、以下のカスタムライブラリへのパスがプログラム内で指定されています。これらは提供されているコードには含まれていませんが、プログラムの実行には必要です。

c:/Programs/python/lib および d:/Programs/python/lib にある tkfile, tkre, tkutils, tksci, tkcrystal

これらのカスタムライブラリは、VASPの出力ファイル解析、結晶構造操作、ファイルI/O、汎用ユーティリティ関数など、本プログラムの主要な機能を提供しています。具体的には、tkVASP はVASP関連ファイルの読み込み、tkCrystal は結晶構造の表現と操作、tkCIF はCIFファイルの生成、tkAtomType は原子種情報の取得などに用いられています。

必要な入力ファイル

プログラムは、VASPの標準的な出力ファイル群を読み込みます。これらは通常、MD計算が実行されたディレクトリ内に存在します。

VASP入力/出力ファイル群:
- OUTCAR: VASPの主要な出力ファイルで、各MDステップの結晶構造（格子ベクトル、原子位置）、エネルギー、温度、原子にかかる力などが記録されています。MDシミュレーションの軌跡の主な情報源となります。
- INCAR: VASPの入力設定ファイル。SYSTEM タグからサンプル名を取得し、POTIM タグからMDステップの時間刻み dt を取得します。
- POSCAR: VASPの初期構造ファイル。MSD計算の基準位置 (origin='POSCAR') や、AXSFファイル作成時の初期構造情報として使用されます。
- CONTCAR: VASPの最終構造ファイル。プログラムはこれを最終構造のCIFファイル出力に使用します。

これらのファイルは、プログラム実行時に引数で指定する CAR_path (通常は OUTCAR ファイルがあるディレクトリ、または OUTCAR ファイル自体) の中に存在することを期待します。

生成される出力ファイル

プログラムは、実行モードと引数に応じて以下のファイルを生成します。sample_name は INCAR の SYSTEM タグから取得されます。

1. `make_axsf` モード (またはMSD計算時にAXSFファイルが存在しない場合)

md.xsf: (例: OUTCAR と同じディレクトリに生成)
- XCrySDenで可視化可能なアニメーションXSF形式のファイル。各MDステップの結晶構造（格子ベクトル、原子番号、カルテシアン座標、原子にかかる力または速度）が格納されています。
inf.csv: (例: OUTCAR と同じディレクトリに生成)
- MDステップごとの時間、温度、エネルギー (TOTEN, EKIN, EKIN_LAT, ETOTAL) などの情報がCSV形式で記録されます。

2. `make_history` モード (MSD計算)

[sample_name]-initial.cif: (例: INCAR と同じディレクトリに生成)
- POSCAR から読み込まれたMD計算の初期結晶構造をCIF (Crystallographic Information File) 形式で保存したもの。
[sample_name]-final.cif: (例: INCAR と同じディレクトリに生成)
- CONTCAR から読み込まれたMD計算の最終結晶構造をCIF形式で保存したもの。
[sample_name]-average.csv: (例: INCAR と同じディレクトリに生成)
- デバッグ目的のCSVファイル。各MDステップにおける系の重心位置、各原子の累積オフセット、および平均原子位置 (初めの nskip ステップを除いた平均) が記録されます。
[sample_name]-history.csv: (例: INCAR と同じディレクトリに生成)
- MDステップごとの詳細な履歴データがCSV形式で記録されます。
  - step: MDステップ数
  - t(fs): 時刻 (フェムト秒)
  - T(K): 温度 (ケルビン)
  - Etot, EKIN, EKIN_LAT, ETOTAL: 各種エネルギー値
  - a, b, c, alpha, beta, gamma, V: 格子定数 (Åおよび度) と単位胞体積 (Å^3)
  - msd([atom_name])(A^2) または msd([site_label])(A^2): 平均二乗変位 (Å^2)。mode オプションに応じて、原子種ごとの平均値、または個々のサイトごとの値が出力されます。

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

(a) XSFファイルの作成 (MD軌跡の可視化用)

python make_md_history.py CAR_path init nmax_read

CAR_path: VASPの出力ファイル (OUTCAR, INCAR, POSCAR など) が存在するディレクトリのパス、または OUTCAR ファイル自体へのパス。
init: XSFファイル作成モードを指定するキーワード。
nmax_read: (オプション) 読み込むMDステップの最大数。デフォルトは-1 (全て読み込む)。

(b) MSDの計算と履歴ファイルの作成

python make_md_history.py CAR_path mode origin drift_correction nskip nmax_read

CAR_path: VASPの出力ファイルが存在するディレクトリのパス、または OUTCAR ファイル自体へのパス。
mode: MSDの計算方法を指定します。
- 'atom' : 原子種ごとにMSDを計算します。
- 'site' : 各原子サイトごとにMSDを計算します。
origin: MSD計算の基準となる初期位置を指定します。
- 'POSCAR' : POSCAR ファイルの初期構造を基準とします。
- 'OUTCAR' : OUTCAR ファイルの最初のMDステップ構造を基準とします。
- 'average' : シミュレーションの最初の nskip ステップを除いた平均位置を基準とします。
drift_correction: 重心ドリフト補正の有無を指定します。
- 'no' : 補正なし。
- 'cm' : 重心移動を補正します。
nskip: (オプション) MSD計算の基準となる位置を平均する際にスキップするMDステップ数。デフォルトは0。
nmax_read: (オプション) 読み込むMDステップの最大数。デフォルトは-1 (全て読み込む)。

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

ここでは、./vasp_md_run というディレクトリにVASPのMD計算結果ファイル (OUTCAR, INCAR, POSCAR, CONTCAR) が格納されていると仮定します。INCAR ファイルには SYSTEM = LiFeO2 と記述されているものとします。

1. XSFアニメーションファイルの作成

VASPのMD軌跡をXCrySDenで可視化するために、md.xsf ファイルを作成します。この例では、最初の2000ステップまでを読み込みます。

python make_md_history.py ./vasp_md_run init 2000

実行結果: このコマンドを実行すると、./vasp_md_run/ ディレクトリ内に以下のファイルが生成されます。

./vasp_md_run/md.xsf: MDシミュレーションの最初の2000ステップまでの結晶構造アニメーションファイル。
./vasp_md_run/inf.csv: 各ステップの温度やエネルギーが記録されたCSVファイル。

2. MSDの計算 (原子種別、重心補正あり)

MD計算の全ステップ (nmax_read を省略するとデフォルトの-1が適用され全て読み込む) について、Li原子とFe原子、O原子それぞれのMSDを計算します。基準位置は最初の nskip=100 ステップを除いた平均位置とし、重心ドリフト補正 (cm) を適用します。

python make_md_history.py ./vasp_md_run atom average cm 100 -1

実行結果: このコマンドを実行すると、まず md.xsf が存在しない場合は自動的に生成されます。その後、./vasp_md_run/ ディレクトリ内に以下のファイルが生成されます。

./vasp_md_run/LiFeO2-initial.cif: POSCAR から生成された初期構造のCIFファイル。
./vasp_md_run/LiFeO2-final.cif: CONTCAR から生成された最終構造のCIFファイル。
./vasp_md_run/LiFeO2-average.csv: 重心位置や各原子の平均位置などが記録されたデバッグ用CSVファイル。
./vasp_md_run/LiFeO2-history.csv: 各MDステップにおける時間、温度、格子情報、そしてLi, Fe, Oそれぞれの原子種のMSDが記録されたCSVファイル。

LiFeO2-history.csv の出力例 (一部):

step,t(fs),T(K),Etot,EKIN,EKIN_LAT,ETOTAL,a,b,c,alpha,beta,gamma,V,msd(Li)(A^2),msd(Fe)(A^2),msd(O)(A^2)
0,0.0,300.0, -1000.0,100.0,50.0,-950.0,5.000,5.000,5.000,90.0,90.0,90.0,125.0,0.000000,0.000000,0.000000
1,1.0,305.2, -999.5,100.5,50.2,-949.3,5.001,4.999,5.001,90.01,89.99,90.00,125.05,0.000001,0.000000,0.000000
...

MSDグラフ作成の考慮事項: 生成された LiFeO2-history.csv の t(fs) 列と msd(AtomName)(A^2) 列をプロットすることで、MDシミュレーション中の各原子種の拡散挙動をグラフで視覚的に分析できます。長時間スケールでMSDが時間に比例して増加する場合、その原子種は拡散していることを示し、その傾きから拡散係数を求めることができます。