CP2KによるMD計算

分子動力学(molecular dynamics: MD)とは

• MDは、原子や分子の動的挙動を時系列で解析する手法。シミュレーションにより、系の構造変化、反応経路、熱的性質などを詳細に解析できる。実験では捉えにくい微視的な現象を数値計算で再現することで、物質の理解を深めるための強力な手段となる。
• MDシミュレーションでは、各時間ステップでエネルギーと力(force)を必要とする。この計算方法により以下に大別できる。
  1. 古典力場を使ってエネルギーと力を計算 -> 古典MD
  2. 量子力学を使ってエネルギーと力を計算 -> 第一原理MD (ab initio MD)
• このチュートリアルは、CP2K を用いた分子動力学シミュレーションの手順および設定例を示す.

サンプルインプットファイル

• 以下はH2O分子に対するMD計算用のサンプルインプットファイルである。各セクションの役割とパラメータを詳細に解説する.

&GLOBAL
  PROJECT water_md      # シミュレーションプロジェクト名を指定
  RUN_TYPE MD          # 実行タイプ。MD は分子動力学を意味する
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE NVT       # NVT アンサンブルを設定（粒子数、体積、温度一定）
    STEPS 5000         # シミュレーションのタイムステップ数
    TIMESTEP 0.5       # 各タイムステップの間隔
    TEMPERATURE 300    # シミュレーション温度（K）
    &THERMOSTAT
      TYPE NOSE        # 温度制御手法としてNOSEサーモスタットを使用
      TIMECON 100      # 温度制御の緩和時間
    &END THERMOSTAT
  &END MD
&END MOTION

&FORCE_EVAL
  METHOD QS           # 力評価法の指定。QSは量子化学計算コード(QuickStep)を表す
  &DFT
    BASIS_SET_FILE_NAME BASIS_SET          # 基底関数セットファイルの名称
    POTENTIAL_FILE_NAME POTENTIAL          # ポテンシャルファイルの名称
    &MGRID
      CUTOFF 400
    &END MGRID
    &XC
      &XC_FUNCTIONAL PBE      # 交換相関（XC）の関数。PBE を指定
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC 10.0 10.0 10.0              # 計算セルの各辺の長さ
      ALPHA_BETA_GAMMA 90.0 90.0 90.0 # セルの角度
    &END CELL
    &COORD
      H  0.0  0.0  0.0              # 原子種とその座標（x, y, z）
      O  0.0  0.0  1.0
      H  1.0  0.0  0.0
    &END COORD
    &KIND H
      BASIS_SET DZVP-MOLOPT-SR-GTH         # 水素用基底関数
      POTENTIAL GTH-POTENTIALS             # 水素用ポテンシャル
    &END KIND
    &KIND O
      BASIS_SET DZVP-MOLOPT-SR-GTH         # 酸素用基底関数
      POTENTIAL GTH-POTENTIALS             # 酸素用ポテンシャル
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

各タグの詳細解説

• &GLOBAL ... &END GLOBAL
  全体設定を記述。

  • PROJECT: プロジェクト名。シミュレーション結果の識別に使用する。
  • RUN_TYPE: 実行するシミュレーションの種類（例: MD）。
  • PRINT_LEVEL: 出力情報の詳細度を設定する.

• &MOTION ... &END MOTION
  運動に関する設定を記述。

  • &MD ... &END MD: 分子動力学固有のパラメータ（エンサンブル、タイムステップ、温度など）を設定.
  • ENSEMBLE: 使用する熱力学アンサンブル（ここでは NVT）を指定.
  • STEPS: シミュレーションの実行ステップ数.
  • TIMESTEP: 各ステップの時間間隔.
  • TEMPERATURE: 目標温度.
  • &THERMOSTAT ... &END THERMOSTAT: 温度制御パラメータを記述.
    • TYPE: 使用するサーモスタットの種類（ここでは NOSE）.
    • TIMECON: 温度制御の緩和時間.

• &FORCE_EVAL ... &END FORCE_EVAL
  力評価とエネルギー計算に関する設定.

  • METHOD: 力計算方法。ここでは QS を指定し、量子力学的な評価を行う.

  • &DFT ... &END DFT: 密度汎関数理論（DFT）による計算設定.

    • BASIS_SET_FILE_NAME と POTENTIAL_FILE_NAME: 使用する基底関数およびポテンシャルのファイル名.
    • &MGRID ... &END MGRID: 密度計算用グリッドのパラメータ.
      • CUTOFF: 計算に使うカットオフ値.
    • &XC ... &END XC: 交換相関計算に関する設定.
      • &XC_FUNCTIONAL ... &END XC_FUNCTIONAL: 具体的な交換相関汎関数（例: PBE）を指定.

  • &SUBSYS ... &END SUBSYS: システム（サブシステム）の構造情報.

    • &CELL ... &END CELL: 計算セルの大きさおよび形状.
      • ABC: セルの各辺の長さ.
      • ALPHA_BETA_GAMMA: セルの角度.
    • &COORD ... &END COORD: 原子の種類と座標情報.
    • &KIND ... &END KIND: 原子種ごとの基底関数およびポテンシャル設定.

実行手順

• 以下のコマンドで計算を実行

cp2k.popt -i cp2k_input.inp -o cp2k_output.out

• cp2k_input.inp: 上記サンプルファイルを保存したファイル名
• cp2k_output.out: シミュレーション出力が保存されるファイル名

結果の解析

• .xyzファイルが出力されるので、可視化する
  • ASEを使うのが楽かもしれない。pip install ase -> ase gui water_md.xyz