make_md_history.py 技術ドキュメント

プログラムの動作

make_md_history.py は、VASP (Vienna Ab initio Simulation Package) の分子動力学 (MD) シミュレーションの出力ファイルを解析し、その結果を標準的なデータ形式（CIFおよびCSV）に変換するPythonスクリプトです。このプログラムの主な目的は、VASPの生出力ファイルから重要な物理量と構造の時系列データを抽出し、他の解析ツールや可視化ツールで利用しやすい形式で提供することです。

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

  • VASP出力ファイルの読み込み: MDシミュレーションのOUTCAR、INCAR、POSCAR、CONTCARファイルを読み込み、シミュレーション情報を抽出します。

  • 構造情報の出力: シミュレーションの初期構造と最終構造をCIF (Crystallographic Information File) 形式で出力します。

  • 時系列データの抽出: 各MDステップにおける時間、温度、全エネルギー、運動エネルギー、格子定数、セル角度、セル体積などの物理量を抽出します。

  • 平均二乗変位 (MSD) 計算: 各原子タイプについて、平均二乗変位 (MSD: Mean Square Displacement) または統合平均二乗変位 (IMSD: Integrated Mean Square Displacement) を計算します。

  • CSV形式での履歴データ出力: 抽出した物理量と計算したMSD/IMSDをCSV (Comma-Separated Values) 形式で出力し、MDシミュレーションの履歴を一元的に記録します。

本プログラムは、VASPのMDシミュレーション結果を手動で解析する手間を省き、結果の可視化や詳細な分析を効率化することを課題解決としています。

原理

本プログラムは、VASPの分子動力学シミュレーションのデータから以下の物理量を抽出し、計算します。

  1. セル情報: 各ステップにおける格子定数 (\(a, b, c\))、セル角度 (\(\alpha, \beta, \gamma\))、およびセル体積 (\(V\))。これらはVASPのOUTCARファイルから直接読み取られます。

  2. 熱力学的量: 各ステップにおける温度 (\(T\))、全エネルギー (TOTEN)、運動エネルギー (EKIN)、格子運動エネルギー (EKIN_LAT)、および総エネルギー (ETOTAL)。これらもOUTCARファイルから抽出されます。

  3. 平均二乗変位 (MSD): MDシミュレーションにおける原子の拡散挙動を評価するための指標です。ある原子タイプに属する \(N\) 個の原子について、時刻 \(t\) における平均二乗変位は以下の数式で定義されます。

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

    ここで、\(\mathbf{r}_i(t)\) は時刻 \(t\) における原子 \(i\) の位置ベクトル、\(\mathbf{r}_i(0)\) は初期時刻 \(t=0\) における原子 \(i\) の位置ベクトルです。プログラムでは、VASPの周期境界条件を考慮して、2つの位置間の最短距離を計算する cry.GetNearestInterAtomicDistance メソッドを用いて \(|\mathbf{r}_i(t) - \mathbf{r}_i(0)|\) を求めています。

  4. 統合平均二乗変位 (IMSD): 本プログラムにおけるIMSDは、各ステップでの「直前のステップからの変位の二乗」を累積したものです。ある原子タイプに属する \(N\) 個の原子について、時刻 \(t\) におけるIMSDは以下のように定義されます。

    \[ \text{IMSD}(t) = \text{IMSD}(t-\Delta t) + \frac{1}{N} \sum_{i=1}^{N} |\mathbf{r}_i(t) - \mathbf{r}_i(t-\Delta t)|^2 \]

    ここで、\(\Delta t\) はMDの時間刻みです。IMSDは、直前のステップからの微小な動きの蓄積を示す指標として機能します。これは、MSDが初期位置からの累積変位を見るのに対し、IMSDは各瞬間の変位を積み上げていく点で異なります。

これらの計算された物理量は、VASPのINCARファイルから読み込まれる時間刻みPOTIMと組み合わせて、時間軸を持つ時系列データとして出力されます。

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

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

  • NumPy: 数値計算を効率的に行うためのPythonライブラリです。

    • インストール方法:

      pip install numpy

  • tklib: プログラム内でtklib.tkfile, tklib.tkutils, tklib.tksci, tklib.tkmatrix, tklib.tkcrystal, tklib.tkapplicationなどのモジュールがインポートされています。

    • このtklibは標準的な公開ライブラリではないため、pipでの直接インストールはできません。このライブラリは、本プログラムの作者によって作成された、VASPの出力解析や結晶構造操作などの機能を提供するカスタムライブラリであると推測されます。

    • インストール方法: tklibは公開されている標準ライブラリではないため、別途入手・インストールが必要です。具体的なインストール方法はtklibの提供元にご確認ください。通常は、tklibのソースコードを適切なPythonのパッケージ検索パスに配置するか、本プログラムと同じディレクトリに配置する必要があります。

必要な入力ファイル

プログラムを実行するには、VASPの分子動力学シミュレーションが完了したディレクトリと、その中に以下のファイルが存在する必要があります。

  1. OUTCAR:

    • VASPの計算結果が出力されるメインファイルです。

    • MDシミュレーションの各ステップにおける原子座標、セル情報、温度、各種エネルギーなどの詳細な情報が記述されています。プログラムはこのファイルからMDの履歴データ（crystal structures セクションなど）を読み取ります。

  2. INCAR:

    • VASPの入力設定ファイルです。

    • プログラムはここからSYSTEMタグ（サンプル名を特定するため）とPOTIMタグ（MDの時間刻みを特定するため）の値を読み取ります。

  3. POSCAR:

    • VASPの初期構造ファイルです。

    • プログラムはシミュレーション開始時の結晶構造を読み取るために使用します。

  4. CONTCAR:

    • VASPの最終構造ファイルです。

    • プログラムはシミュレーション終了時の結晶構造を読み取るために使用します。

これらのファイルは、指定されたoutcar_path引数（VASP出力ディレクトリまたはOUTCARファイルへのパス）から見つけ出されます。

生成される出力ファイル

プログラムの実行により、以下のファイルが生成されます。

  1. [sample_name]-initial.cif:

    • ファイル名: INCARファイル内のSYSTEMタグの値（例: MgO）がsample_nameとして使用されます。

    • 内容: VASPのPOSCARから読み取られた、シミュレーション開始時点の結晶構造情報がCIF (Crystallographic Information File) 形式で保存されます。

  2. [sample_name]-final.cif:

    • ファイル名: INCARファイル内のSYSTEMタグの値（例: MgO）がsample_nameとして使用されます。

    • 内容: VASPのCONTCARから読み取られた、シミュレーション終了時点の結晶構造情報がCIF形式で保存されます。

  3. [sample_name]-history.csv:

    • ファイル名: INCARファイル内のSYSTEMタグの値（例: MgO）がsample_nameとして使用されます。

    • 内容: MDシミュレーションの時系列履歴データがCSV形式で保存されます。各行は1つのMDステップに対応し、以下の情報が含まれます。

      • step: MDステップの数。

      • t(fs): シミュレーション経過時間 (フェムト秒)。

      • T(K): 温度 (ケルビン)。

      • Etot: 全エネルギー。

      • EKIN: 運動エネルギー。

      • EKIN_LAT: 格子運動エネルギー。

      • ETOTAL: 総エネルギー。

      • a,b,c: 格子定数 (Å)。

      • alpha,beta,gamma: セル角度 (度)。

      • V: セル体積 (Å^3)。

      • [msd/imsd](AtomType1)(A^2): 各原子タイプごとの平均二乗変位 (または統合平均二乗変位) (\(Å^2\))。mode引数によってmsdまたはimsdのどちらかになります。

  4. make_md_hisotry-out.txt:

    • ファイル名: このプログラム自体の実行ログファイルです。

    • 内容: プログラムの実行中に標準出力に表示される全ての情報（読み込みファイルパス、計算状況、エラーメッセージなど）がここにリダイレクトされて記録されます。このファイルは、outcar_pathで指定されたディレクトリ内に作成されます。

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

make_md_history.py は、以下の基本構文でコマンドラインから実行します。

python make_md_history.py [mode] [outcar_path]

  • [mode] (必須): 計算する平均二乗変位の種類を指定します。

    • msd: 標準の平均二乗変位 (Mean Square Displacement)。各ステップにおいて、初期位置からの変位の二乗を原子タイプごとに平均します。

    • imsd: 統合平均二乗変位 (Integrated Mean Square Displacement)。各ステップにおいて、直前のステップからの変位の二乗を原子タイプごとに平均し、それを累積していきます。

  • [outcar_path] (必須): VASPの出力ファイル群（OUTCAR, INCAR, POSCAR, CONTCAR）を含むディレクトリのパス、またはOUTCARファイル自体へのパスを指定します。プログラムは指定されたパスからこれらのファイルを自動的に探索します。

例:

python make_md_history.py msd ./path/to/vasp_md_run

この例では、./path/to/vasp_md_run ディレクトリ内のVASP出力ファイルを解析し、MSDを計算してCSVファイルに出力します。

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

ここでは、make_md_history.py を使用する具体的なシナリオと、その実行結果の例を示します。

シナリオ: VASPのMDシミュレーションをディレクトリ vasp_md_data で実行したとします。このディレクトリには、OUTCAR, INCAR, POSCAR, CONTCAR が含まれており、INCAR ファイルには以下の設定が含まれているとします。

SYSTEM = "Example-System"
POTIM  = 1.5

また、POSCARにはMgとOの2種類の原子が含まれているとします。

実行コマンド: ここでは、直前ステップからの変位を累積する統合平均二乗変位 (imsd) を計算します。

python make_md_history.py imsd vasp_md_data

実行結果の説明:

  1. 標準出力 (および vasp_md_data/make_md_hisotry-out.txt へのリダイレクト): プログラムは、解析対象のディレクトリ、読み込むファイル、抽出されたSYSTEM名やPOTIM値、そして各MDステップでのセル情報（格子定数、角度）などを逐次出力します。

    Open logfile [vasp_md_data/make_md_hisotry-out.txt]

=============== Convert VASP output files to XCrySDen AXSF file ============

mode   :  imsd
CAR dir:  vasp_md_data
INCAR  :  vasp_md_data/INCAR
POSCAR :  vasp_md_data/POSCAR
CONTCAR:  vasp_md_data/CONTCAR
OUTCAR :  vasp_md_data/OUTCAR

*** Read initial crystal structure from [vasp_md_data/POSCAR]
cell:   3.10000000   3.10000000   3.10000000 A     90.000000   90.000000   90.000000

*** Read final crystal structure from [vasp_md_data/CONTCAR]
cell:   3.10500000   3.10500000   3.10500000 A     90.000000   90.000000   90.000000

Read [vasp_md_data/INCAR]

sample_name: Example-System
POTIM: 1.5 fs

Write to [Example-System-initial.cif]
Write to [Example-System-final.cif]

Write history to [Example-System-history.csv]
# of crystal structures: 1000

# of atom sites for types =  [4 4]  (例: Mg原子が4つ、O原子が4つ)

Read crystal structures from [vasp_md_data/OUTCAR]
step #00000:   3.10000000   3.10000000   3.10000000 A     90.000000   90.000000   90.000000
step #00001:   3.10012345   3.10012345   3.10012345 A     90.000010   90.000010   90.000010
...

  2. 生成されるファイル: vasp_md_data ディレクトリ内に以下の3つのファイルが生成されます。

    • Example-System-initial.cif: MD開始時の結晶構造。

    • Example-System-final.cif: MD終了時の結晶構造。

    • Example-System-history.csv: MD履歴データ。

  3. Example-System-history.csv の内容例:

    step,t(fs),T(K),Etot, EKIN, EKIN_LAT, ETOTAL,a,b,c,alpha,beta,gamma,V,imsd(Mg)(A^2),imsd(O)(A^2)
0,0.0,300.00,  -10.1234,   1.2345,   1.2345,  -8.8889,3.10000000,3.10000000,3.10000000,90.000000,90.000000,90.000000,29.79100000,0.00000000,0.00000000
1,1.5,301.25,  -10.1230,   1.2350,   1.2350,  -8.8880,3.10012345,3.10012345,3.10012345,90.000010,90.000010,90.000010,29.80012345,0.00123456,0.00078901
2,3.0,300.98,  -10.1232,   1.2348,   1.2348,  -8.8884,3.10024690,3.10024690,3.10024690,90.000020,90.000020,90.000020,29.80924690,0.00250000,0.00160000
...