水和構造を作成し分子動力学シミュレーションを実行してみる

# 概要
本稿ではタイトルの通り水和構造を作成し，分子動力学シミュレーションを実行する手続きを紹介します．水和構造はpackmol プログラムを用いて作成します．分子動力学シミュレーションはLAMMPSコードを用います．

# 用いるアプリケーション
本稿の作業においては様々なアプリケーションを利用して目的の処理を行います．それぞれについて簡単に説明したいと思います．

- 水和構造の作成：[packmol](http://leandro.iqm.unicamp.br/m3g/packmol/home.shtml); 分子を指定の領域以内に指定の数だけランダムに配置することのできるアプリケーションです．
- 座標データ変換：[conv.py](https://azuma.nims.go.jp/); いろいろな形式の座標データファイルを変換することができるツールです．第一原理擬ポテンシャルプログラム[PHASE/0](https://azuma.nims.go.jp/software/phase)に同梱されています．
- 原子座標データの可視化：[Jmol](http://jmol.sourceforge.net/); Javaで動作する原子配置可視化アプリケーションです．動作が軽快です（個人の感想）
- 分子動力学シミュレーター : [LAMMPS](https://www.lammps.org/); 汎用の分子動力学シミュレーターです．多くの原子間ポテンシャルをサポートします．
- mp4ファイルの作成 : [OpenCV](https://opencv.org/); 画像や動画データをプログラム的に扱う際に便利な機能を提供してくれるライブラリーです．
- animated GIFの作成：[ImageMagick](https://imagemagick.org/index.php); 画像変換などを行うことができるアプリケーションです．animated GIFの作成にはImageMagickのconvertコマンドを使います．

# 水和構造の作成
水和構造の作成にはpackmolプログラムを利用します．系は何でもよいのですが，二酸化酸素分子が水に溶けている問題を採用したいと思います．
## packmolの入手・コンパイル
packmolをコンパイルするにはFortranコンパイラーが必要です．お使いのシステムにインストールされていない場合インストールしておく必要があります．

packmolは配布元のウェブサイトから入手することができます．本稿執筆時の最新版，バージョン20.3.5を用いることにします．zipもしくはtar.gzファイルをダウンロードし，好きな場所に保存しましょう．解凍後以下の手続きでコンパイルすることができます．
```
$ cd packmol-20.3.5
$ ./configure
$ make
```

## packmolのインプットの作成
ここでは三辺12Åの立方体に水分子60個，二酸化炭素分子10個が詰められた構造を作成します．そのためのpackmolの入力は以下のようになります．
```
tolerance 2.0
filetype xyz
output out.xyz
 
structure H2O.xyz
  number 60
  inside box 0 0 0 10 10 10
end structure

structure CO2.xyz
  number 10
  inside box 0 0 0 10 10 10
end structure
```
`tolerrence 2.0`とすると分子間の距離は最低でも2Åは確保されることになります．`filetype xyz`は採用する座標データファイル形式を指定するものであり，xyzというのは後で紹介する形式です．`output out.xyz`によって結果を`out.xyz`というxyzファイルに出力することが指定されています．`structure ... end structure`によって水和構造の種となる分子の情報を指定しています．`structure`に続く文字列によって分子座標が記録されたファイルを指定します．`number`によってその分子の水和構造における個数を指定します．`inside box 0 0 0 10 10 10`によって分子を三辺10Åの立方体の中に詰めることを指定します．最終的な単位胞は三辺12Åの立方体ですが，packmolは周期境界条件を考慮した上で分子を詰めることはしてくれず，分子が重ならないようにするためには少し余裕をもたせる必要があります．この余裕分は，数千ステップ分の分子動力学シミュレーションによって解消することができます．
## 分子座標データ
ここでは分子座標を指定するファイルとして`H2O.xyz`と`CO2.xyz`というファイルを用います．前者が水分子，後者が二酸化炭素分子に対応します．その内容は以下の通り．

### H2O.xyz
```
3
# H2O
H        0.95   -0.55   0
H       -0.95   -0.55   0
O       0       0       0
```

### CO2.xyz
```
3
# CO2
O        1.4    0       0
O       -1.4    0       0
C       0       0       0
```

## packmolの実行
以下の要領でpackmolを実行します（ホームディレクトリーにpackmolを展開したと仮定しています）
```
$ $HOME/packmol-20.3.5/packmol < packmol.in
``` 
しばらくして以下のようなログが得られれば無事水和構造が作成できています．
```
################################################################################

                                 Success!
              Final objective function value: .46489E-01
              Maximum violation of target distance:   0.000000
              Maximum violation of the constraints: .82269E-02

--------------------------------------------------------------------------------

              Please cite this work if Packmol was useful:

           L. Martinez, R. Andrade, E. G. Birgin, J. M. Martinez,
         PACKMOL: A package for building initial configurations for
                   molecular dynamics simulations.
          Journal of Computational Chemistry, 30:2157-2164,2009.

################################################################################
```
得られた構造は`out.xyz`ファイルにxyz形式で記録されます．out.xyzをJmolを用いて可視化すると以下のようになります．
![out.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/262037/f7482b23-84c9-d63b-cb29-a36b822ef7c2.jpeg)

# 分子動力学シミュレーションの実行
つぎに，作成した水和構造を初期配置とした分子動力学シミュレーションを実行してみましょう．用いる分子動力学シミュレーターはLAMMPS, 原子間ポテンシャルとしては（セットアップが易しい）ReaxFFを用いることにします．

## 座標データの変換
作成した水和構造はxyz形式で出力されていますが，これをLAMMPSに直接読み込ませることはできません．そこで，LAMMPSが読み込める形式に変換する作業を行います．いろいろな方法がありますが，ここではPHASE/0プログラムに付属する座標データコンバーター`conv.py`を用いることにします．

xyzファイルには単位胞の情報は記録されませんが，古典分子動力学シミュレーションを行うためには単位胞の情報を与える必要があります．そこで，変換を行う前に`conv.py`が処理できる形式でxyzファイルに単位胞の情報を記録します．具体的にはxyzファイルの二行目に$a \ b \ c \  \alpha \ \beta \ \gamma$という形式で格子定数を指定すれば`conv.py`は単位胞の情報を含む座標データファイルを作成してくれます．長さの単位はÅ, 角度の単位は度単位です． 適当なテキストエディターを用いて編集するか，以下のコマンドを利用してください．
```
$ sed --in-place "2c 12 12 12 90 90 90" out.xyz
```
格子定数の書き込みができたら，`conv.py`を以下のように実行します（PHASE/0プログラムがホームディレクトリーにインストールされていることを仮定しています）
```
$ $HOME/phase0_2021.02/bin/conv.py --batch --input_type=xyz --output_type=LAMMPS_input --input_file=out.xyz --output_file=input.lmp
```
`conv.py`に渡しているオプションについて説明します．``--batch``は「バッチモード」で利用することを意味します．この指定がないと対話的に利用することになります．``--input_type=xyz``によって変換元のファイルの形式がxyz形式であることを指定しています．``--output_type=LAMMPS_input``によって変換先のファイルの形式がLAMMPSの入力形式であることを指定します．``--input_file=out.xyz``で変換元のファイル名を，``--output_file=input.lmp``によって変換先のファイル名を指定しています．
## 電荷に相当するカラムの追加
`conv.py`が出力するLAMMPSの入力座標データは最も単純な [atom_style atomic](https://docs.lammps.org/atom_style.html)に対応したものです．今回用いる予定のReaxFFポテンシャルはこれよりももう少しややこしく，各原子の電荷の指定が必要な[atom_style charge](https://docs.lammps.org/atom_style.html)なので，編集が必要です．

`atomic`スタイルは`ID elemenentID x y z`のような形式で原子座標を指定します．これに対しchargeスタイルは`ID elementID charge x y z`のような形式で原子座標を指定します．電荷はReaxFFの場合平衡化によってほどよい値を算出してくれるので初期値は重要ではありませんが，入力のカラムとしては定義されている必要があります．ここではとりあえず0に設定することにします．以下の要領で編集することができます．
```
$ sed "s/ 1 / 1 0 /g" input.lmlp | sed "s/ 2 / 2 0 /g" | sed "s/ 3 / 3 0 /g" > input_q.lmp
```
`スペース整数値スペース`という並びは元素指定のカラムの前後のみにしか現れないので，sedコマンドによってそれを検索し，`スペース整数値スペース0スペース`という並びに置換しています．この置換処理を3種類ある元素すべてに対して行い，結果を`input_q.lmp`というファイルに記録しています．

## LAMMPSのインプットスクリプト
LAMMPSを用いて分子動力学シミュレーションを行うためには，座標データファイルのほかインプットスクリプトを用意する必要があります．今回は以下のようなものを利用することにします．
```
units         real
boundary p p p
atom_style    charge
read_data     input_q.lmp
mass 1 1
mass 2 16
mass 3 12
pair_style    reax/c NULL
pair_coeff    * * ffield.reax.cho H O C
timestep      0.25
dump          1 all custom 200 dump element type xsu ysu zsu
dump_modify   1 element H O C
thermo 100
fix           2 all qeq/reax 1 0.0 10.0 1e-6 reax/c
minimize        1.0e-6  1.0e-8  1000  10000
velocity all  create 300.0 6345
unfix 2
fix           1 all nve
fix           2 all qeq/reax 1 0.0 10.0 1e-6 reax/c
fix           3 all temp/berendsen 300.0 300.0 100.0
run           10000
```
`units real`で単位系を指定しています．`boundary p p p`によって三方向周期的であることを指定しています．`atom_style charge`で用いる`atom_style`の指定を行っています．`read_data input_q.lmp`によって先ほど作成した座標データファイルを読み込む指示をしています．`mass 1 1`などの行で各元素の質量を指定しています．`pair_style reax/c NULL`によってReaxFFポテンシャル使用を指定し，さらに`pair_coeff    * * ffield.reax.cho H O C`によってReaxFFの設定を行っています．ReaxFFを用いるには対応するポテンシャルファイルが必要ですが，この例ではLAMMPSに同梱されている`ffield.reax.cho`というファイルを用いています．また，元の原子の並びから元素1が水素，元素2が酸素，元素3が炭素になっていますが，これをReaxFFのポテンシャルファイルの内容とマップするために`ffield.reax.cho H O C`としています．`timestep 0.25`で時間きざみを0.25 fsとしています．`dump          1 all custom 200 dump element type xsu ysu zsu`でdumpファイル（LAMMPSにおいて座標データ履歴を記録するファイル）の設定を行っています．後で`conv.py`を用いて変換することを見越した設定が施されています．この行の200という数値は200MDステップに一度座標データをdumpファイルに記録することを意味します．また，dumpという文字列はdumpファイルのファイル名としてdumpが用いられることを意味します．`dump_modify   1 element H O C`によってdumpファイルにおいて利用する原子種を表す文字列として単なる整数ではなく元素名が採用されるようになります．`thermo 100`によって温度やエネルギーなどの標準出力への出力頻度を指定しています．このように指定すると，100 MDステップに一度ログに情報が記録されます．`fix           2 all qeq/reax 1 0.0 10.0 1e-6 reax/c`はReaxFFの電荷平衡化の設定です．分子動力学シミュレーションに入る前に構造最適化を行うコマンドが`minimize        1.0e-6  1.0e-8  1000  10000`です．`velocity all  create 300.0 6345`によって温度300Kに相当する初期速度を各原子に割り当てています．unfix 2でいったん先の電荷平衡化の設定を無効にしたあと本番のMDの設定を行います．`fix           1 all nve`によって基本の統計アンサンブルがNVEアンサンブルであることを指定しています．しかしながら，`fix           3 all temp/berendsen 300.0 300.0 100.0`によってBerendsen熱浴によって実際には温度制御を行うことが指定されています．最後に`run           10000`とすることによって一万ステップの分子動力学シミュレーションを実行します．

## LAMMPSの実行
LAMMPSは，たとえば以下の要領で実行します．
```
$ mpiexec -n  lmp_mpi < input.lmp
```
具体的なコマンドは処理系によって異なりますが，入力スクリプトファイル(`input.lmp`)をリダイレクトによってLAMMPSのバイナリーに食わせる記法は共通です．各ステップにおけるエネルギーや温度などのログは標準出力と`log.lammps`ファイルに記録されます．座標データの履歴は`dump`ファイルに出力されます．

# 分子動力学シミュレーションの結果の解析

## エネルギーおよび温度の履歴
エネルギーや温度などの履歴は`log.lammps`ファイルに記録されます．初期化などが行われ，構造最適化も終了すると以下のように指定のステップ数に一度情報が出力されます．
```
Step Temp E_pair E_mol TotEng Press
     156          300   -19750.033            0   -19563.136   -5175.4269
     200    250.56565   -19715.944            0   -19559.844     12992.05
     300    343.00041   -19774.981            0   -19561.296   -11328.316
     400     347.8871   -19783.916            0   -19567.186   -14531.112
     500    370.69045   -19805.794            0   -19574.858   -11352.457
     600    362.30115   -19812.109            0     -19586.4   -11739.861
     700    368.82392   -19826.301            0   -19596.527   -7999.6108
     ...
     ...
```
1カラム目がMDのステップ数，2カラム目が温度，3カラム目がエネルギー，4カラム目がReaxFFでは使われない分子間エネルギー，5カラム目が全エネルギー，6カラム目が圧力に対応します(1カラム目の最初の数値が156なのはまず構造最適化が行われたため) エネルギーや温度とMDステップの関係をプロットしたのが以下の図です．
![e.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/262037/d7efb7fe-37e0-382f-22b0-34ef3dbad516.png)
おおよそ4000ステップくらいまでは大きな変動をみせていますが，それ以降は落ち着いた振る舞いになっています．これは，packmolで初期構造を作る際に間隙を設けたこともあり最初の構造は不安定だったのが，シミュレーションによって落ち着いてきたことを表しています．

## 座標の履歴のアニメーション
dumpファイルから座標の履歴のアニメーションを作る方法を説明します．二つの方法を紹介しますが，いずれもdumpファイルの各フレームの原子配置を連番付きの画像データファイルに保存し，それを連結するものです．
### 連番付き画像ファイルの作り方
連番付きの画像ファイルを作成する方法はいろいろあると思いますが，ここではJmolを使った方法を紹介します．まずはdumpファイルをJmolが読める形式に変換します．これまでと同様`conv.py`を用いて`xyz`形式に変換していきます．
```
$ $HOME/phase0_2021.02/bin/conv.py  --batch --input_type=LAMMPS_output --output_type=xyz --input_file=dump --output_file=dump.xyz
```
`--input_type=LAMMPS_output`で変換元がLAMMPS_output形式，すなわちdumpファイルであること，`--output_type=xyz`で変換先がxyz形式であることを指定しています．`--input_file=dump`と`--output_file=dump.xyz`によって変換元および変換先のファイルのファイル名がそれぞれ`dump`および`dump.xyz`であることを指定しています．このコマンドによって`dump.xyz`というxyz形式の座標データファイルに変換されます．

ついでJmolで作成したxyz形式ファイルを開きます．これはJmolへの引数として`dump.xyz`を渡すか，「ファイル」メニューから「開く」を選択し，結果得られるファイル選択ダイアログを介して指定してください．構造が読み込まれ，表示されたら「ファイル」メニューから「コンソール」を選びます．結果現れるコンソールから以下のようなコマンドを実行すると作業ディレクトリーに連番付きの画像ファイルが作成されます．このコマンドの`movie`部分は好きな文字列を使ってください．
```
$ write frames {*} 'movie.jpg'
```
このコマンドによって`movie0001.jpg, movie0002.jpg, ...` といった連番付きの画像ファイルが作成されるはずです．

### mp4ファイルの作り方
jpgファイルをつなげてmp4ファイルを作る方法として，[OpenCV](https://opencv.org/)を用いる方法を紹介します．ここではOpenCVのPythonインターフェースを用います．お使いのPythonインタープリターにインストールされていない場合，pipやcondaなどでインストールします．
```
$ pip install opencv-python もしくは
$ conda install opencv
```

連番付き画像データファイルからmp4を作るスクリプト例を以下に記します．
```python
import cv2
import glob

fourcc = cv2.VideoWriter_fourcc('m','p','4','v')
video = cv2.VideoWriter('video.mp4', fourcc, 20, (800, 600))
jpgfiles = glob.glob('*.jpg')
jpgfiles.sort()
for jpgfile in jpgfiles:
    img = cv2.imread(jpgfile)
    img = cv2.resize(img, (options.width,options.height))
    video.write(img)
video.release()
```
OpenCVの機能にアクセスするため，`cv2`というモジュールを`import`しています．`cv2`の`VideoWriter_fourcc`という関数を用いて[fourcc](https://www.fourcc.org/)を作成し，それを`cv2`の`VideoWriter`クラスのコンストラクターに渡しています．VideoWriterにはほかにも結果を出力するファイルのファイル名(video.mp4)，フレームレート(20)，幅と高さ((800,600))なども渡しています．`glob.glob`によって実行ディレクトリー以下のすべてのjpgファイルのファイル名を取得し，ソートしたあとforループをまわしています．`cv2.imread`でjpgファイルを読み込み，さらに`cv2.resize`でリサイズしたあと`video.write`によって一コマずつ動画にフレームを追加しています．最後に`video.release()`を呼び，終了しています．結果は`video.mp4`というmp4ファイルに記録されます．

### animated GIFの作り方
連番付きの画像データファイルからanimated GIFを作る方法として，[ImageMagick](https://imagemagick.org/index.php)プログラムに付属するconvertコマンドを利用する方法を紹介します．以下の要領でconvertコマンドを実行します．
```
$ convert -delay 10 -loop 0 movie*.jpg  movie.gif
```
`-delay`にはコマの間の間隔を指定します．`-loop`には動画のループ回数を指定します．0とするといつまでもループするanimated GIFが作成されます．`movie*.jpg`とすることによって入力の画像データファイルを複数指定しています．最後の引数に結果を出力するファイル形式とファイル名をmovie.gifという文字列によって指定しています．この処理によって作成されるanimated GIFとは以下のようなものです．

![movie.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/262037/c115b880-2684-60b5-d42c-c29d6e4193a9.gif)

# 終わりに
簡単な手続きによって分子の水和構造を作成し，LAMMPSによって分子動力学シミュレーションを行う例を紹介しました．この例で行ったのは常温かつごく短時間のシミュレーションなのでこれといった興味深い現象は現れませんが，より長時間もしくは高温でシミュレーションを行えば何らかの反応（たとえば炭酸ができるなど）をみることもできるかもしれません．同じ手続きで好きな分子の組み合わせの構造を作成することができますし，少し工夫すれば固液界面をモデリングすることもできます．興味があればチャレンジしてみてください．