概要
本稿ではタイトルの通り水和構造を作成し，分子動力学シミュレーションを実行する手続きを紹介します．水和構造はpackmol プログラムを用いて作成します．分子動力学シミュレーションはLAMMPSコードを用います．
用いるアプリケーション
本稿の作業においては様々なアプリケーションを利用して目的の処理を行います．それぞれについて簡単に説明したいと思います．
- 水和構造の作成：packmol; 分子を指定の領域以内に指定の数だけランダムに配置することのできるアプリケーションです．
- 座標データ変換：conv.py; いろいろな形式の座標データファイルを変換することができるツールです．第一原理擬ポテンシャルプログラムPHASE/0に同梱されています．
- 原子座標データの可視化：Jmol; Javaで動作する原子配置可視化アプリケーションです．動作が軽快です（個人の感想）
- 分子動力学シミュレーター : LAMMPS; 汎用の分子動力学シミュレーターです．多くの原子間ポテンシャルをサポートします．
- mp4ファイルの作成 : OpenCV; 画像や動画データをプログラム的に扱う際に便利な機能を提供してくれるライブラリーです．
- animated GIFの作成：ImageMagick; 画像変換などを行うことができるアプリケーションです．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を用いて可視化すると以下のようになります．
分子動力学シミュレーションの実行
つぎに，作成した水和構造を初期配置とした分子動力学シミュレーションを実行してみましょう．用いる分子動力学シミュレーターは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に対応したものです．今回用いる予定のReaxFFポテンシャルはこれよりももう少しややこしく，各原子の電荷の指定が必要なatom_style chargeなので，編集が必要です．
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ステップの関係をプロットしたのが以下の図です．
おおよそ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を用いる方法を紹介します．ここではOpenCVのPythonインターフェースを用います．お使いのPythonインタープリターにインストールされていない場合，pipやcondaなどでインストールします．
$ pip install opencv-python もしくは
$ conda install opencv
連番付き画像データファイルからmp4を作るスクリプト例を以下に記します．
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を作成し，それを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プログラムに付属するconvertコマンドを利用する方法を紹介します．以下の要領でconvertコマンドを実行します．
$ convert -delay 10 -loop 0 movie*.jpg movie.gif
-delayにはコマの間の間隔を指定します．-loopには動画のループ回数を指定します．0とするといつまでもループするanimated GIFが作成されます．movie*.jpgとすることによって入力の画像データファイルを複数指定しています．最後の引数に結果を出力するファイル形式とファイル名をmovie.gifという文字列によって指定しています．この処理によって作成されるanimated GIFとは以下のようなものです．
終わりに
簡単な手続きによって分子の水和構造を作成し，LAMMPSによって分子動力学シミュレーションを行う例を紹介しました．この例で行ったのは常温かつごく短時間のシミュレーションなのでこれといった興味深い現象は現れませんが，より長時間もしくは高温でシミュレーションを行えば何らかの反応（たとえば炭酸ができるなど）をみることもできるかもしれません．同じ手続きで好きな分子の組み合わせの構造を作成することができますし，少し工夫すれば固液界面をモデリングすることもできます．興味があればチャレンジしてみてください．