Molecular Dynamics by EXCEL

Abstract

EXCELによる分子動力学プログラムを開発した．採用した分子模型はレナードージョーンズ関数である．実行方法，結果の読み方，設定条件の入力の仕方などを示した．ワークシートから系の状態を規定する温度と数密度を与える．熱力学量の他にニ体相関関数，積算配位数，平均ニ乗変位，速度自己相関関数，自己拡散係数などが計算できる．プログラムはEXCEL worksheet として付録に示される．

1 はじめに

分子動力学を手近に体験できるようにEXCELによる分子動力学プログラムを開発した．

使用言語はEXCEL-VBAであるので，EXCELのオプションからVBAの使用を許す設定をすれば容易に実行できる．分子動力学シミュレーションについては書籍 [1]を参考にできる．

2 モデル

12–6レナードージョーンズ相互作用 [2]する球形分子を仮定している．球形に近似できれば多くの分子系に適用できる [3]．相互作用の関数形は，分子間距離rを変数として次の形である．   

$u\left(r\right)=4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]$(1)

ここで$\epsilon$はポテンシャルの深さを示すパラメータであり，$\sigma$は分子直径を表すパラメータである．これらの値を指定すれば具体的な分子を特定できる．パラメータの値は広く知られている [2,3]

プログラムはこれらのパラメータを単位としたときの換算変数で書かれている．たとえば温度Tの換算変数T^rはkをボルツマン定数として次のように定義される．   

$T^r\equiv T/(\epsilon /k)$(2)

同様に時間の単位として次のように定義された$\tau$を使用して，速度は換算されている．なおこの式におけるmは粒子の質量である．   

$\tau =\sqrt{\frac{m{\sigma }^2}{\epsilon }}$(3)

  

$v^r=\frac{v}{\sigma /\tau }$(4)

プログラムではこの添え字rは省略されている．出力では極力$T/(\epsilon /k)$の形で成されるが, 単にTの表記がなされることがある．

3 ダウンロードと解凍

プログラムは本論文の付録として添付されているので，J-STAGE [4]から次のようにダウンロードする．本論文を選択した後，電子付録のタブでDownloadアイコンをクリックする．

Zipファイルがダウンロードできたら，これをダブルクリックで解凍する．解凍できたファイルは自身の作業用フォルダーへコピーする．

4 EXCELの設定

ファイルは固体・液体・気体の典型的なシミュレーションを直ちに実行できるように3つ用意されているので,たとえばNTVsolidN = 256.xlsmを開く．EXCELのタブに「開発」が見当たらない場合は，EXCELの「オプション」から「リボンのユーザ設定」で「メインタブ」の中で「開発」にチェックを入れる．つぎに開発タブを開き，「マクロのセキュリティ」を選択する．「VBA プロジェクト オブジェクト モデルへのアクセスを信頼する」にチェックを入れる．

5 実行の方法

sheet1に入力データの説明がある．"設定条件の数値をワークシートC列のC5, C6・・・に入力してください""入力データの編集が済んだらマクロ「MainMD」を実行してください．"とある．計算を実行するとファイルは上書きされるので，別名を付けて保存を行った後実行する．初めての練習では付録のファイルにある設定条件は変更しないで，開発タブから「マクロ」を開き，「MainMD」を選択し，実行ボタンを押す．N = 256のファイルでは，基本セル内の分子の個数は256で所要時間は12分から5分程度である．計算途中に「応答なし」の表示が出ても, 12分は待つべきである．正常終了の場合はsheet4に「正常終了」と表示される．

6 結果の読み方

多くの結果はsheet1に表示される．sheet2には2体相関関数・積算配位数が示される．

各ステップでの運動エネルギー ek/ε，ポテンシャルエネルギー ep/ε，ハミルトニアンh/ε，ビリアル項vir/ε，圧力 p(ε/σ³)，温度 T(ε/k)などとともに周期境界条件のもとでの粒子1の座標 X (1)/σ，基本セル [1]に戻さない座標 Xr (1)/σ，速度 Vx (1)/(σ/τ)などが1000ステップ分について表として示されている．Figure 1とFigure 2に例を挙げた．

Figure 1.

 The coordinate in the periodic boundary condition X (1)/σ vs. step.

Figure 2.

 The coordinate Xr (1)/σ vs. step.

Figure 1に示した座標は，基本セルからはみ出したとき，周期境界条件 [1]のもとで基本セル内に戻されている様子を示している．Figure 2では平均2乗変位 [1]の計算に使用するために，粒子を基本セル内に戻していない．

7 ニ体相関関数

sheet2にはニ体相関関数 [1] gr (r)と積算配位数 [1]rcn (r)が示されている．ニ体相関関数は通常g (r)と書かれる．本プログラムの中ではgr (r)という名前で書かれている．この例では粒子数N = 256, 数密度 = 1/${\sigma }^3$なので，体積V = 256${\sigma }^3$,基本セルの一辺の長さL = 256^1/3$\sigma$ = 6.35 $\sigma$となる．Figure 3に例を示した．

Figure 3.

 The pair correlation function gr (r) and running coordination number rcn (r) vs. molecular distance r/σ.

本プログラムでは力の計算でminimum image convention [5] を使用している．これは周期境界条件のもとでの相互作用エネルギーと力の計算において，中心に考える分子から見て最も近い距離にある相手分子を一つだけ考慮するものである．分子間距離の統計もこの方法を使用しているため，分子間距離がL/2を超えたときのニ体相関関数と積算配位数は通常の意味と異なる．そこで$r\le L/2$の部分のみを使用する．$r>L/2$の部分も示されているのはプログラムのチェックのためである．

8 平均ニ乗変位

sheet1の下の方には平均ニ乗変位msd(t)と速度自己相関関数c (t)が示されている．平均ニ乗変位msd (t)の例をFigure 4に示した．またFigure 5に速度自己相関関数c (t)と積算数s (t)も示した．積算数とは平均ニ乗変位と速度自己相関関数の計算でとった積算数である．時間tが増えると積算数は減少するので，出力には長時間にわたる関数値が示されているが，時間領域については初めの25%程度までが信頼できると見られる．

Figure 4.

 The mean square displacement msd (t)/σ² vs. time t/τ.

Figure 5.

 The velocity autocorrelation function c (t)/(σ²/τ²) and accumulation number s (t) vs. time t/τ.

9 計算条件の入力法

sheet1の左上の隅(C列)で計算条件を入力する．温度Tempと NumberDensityN/V で分子系の状態(気体・液体・固体)を指定する．つまりレナードージョーンズ系について知られている相図 [6]を参考にする．基本セルに含まれる分子数Nはnfcc で指定する．N = nsp = 4*nfcc³ の関係がある．Nは最低でも256が望ましい．できれば500や864の方が望ましい．1000ステップの計算にかかった時間をsheet4に示してある．ステップ数nstepは最低で1000が必要となる．これを増やすと精度が増すが，計算時間と必要メモリが増大するのでプログラム内で最高で10000と設定してある．運動を数値的に解く刻み幅は0.01のままでよい．初期分子配置における最大歪み幅 deltaは気体で0.3，液体で0.1，固体で0.01程度が望ましい．初期配置は面心立方格子点から，座標を最大歪み幅に範囲[-1, 1]の一様算数を掛けた値だけずらせている．

10 プログラム

プログラムのコードは開発からVisual Basicを選択すると見ることができる．EXCEL VBAについての書籍も多数出版されている．学ぶ場合は入門書 [7]が望ましい．プログラムの構造については別途解説する [8]．プログラムは片岡の著書に基づいている [9]．

References

• [1]   I. Okada, E. Osawa, (ed),”Bunshi Simyureshon Nyumon”, Kaibundo (1989)
• [2]   P. W. Atkins, Physical Chemistry, Oxford Univ. Press, Oxford (1998).
• [3]   F. Cuadros, I. Cachadiña, W. Ahamuda, Mol. Eng., 6, 319 (1996).
• [4]   https://www.jstage.jst.go.jp/browse/jccj/-char/ja/
• [5]   M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford (1992).
• [6]   Y. Kataoka, Y. Yamada, J. Comput. Chem. Jpn., 13, 130 (2014).
• [7]   K. Terasaka,”Kagakukei Gakuseinotameno Excel/VBA Nyumon”, Koronasha (2009).
• [8]   Y. Kataoka, Y. Yamada, Bull. Res. Cent. Comp. Multimedia Stud. Hosei Univ., (2015), (to be published).
• [9]   Y. Kataoka,”Bunshidorikigaku To Monte Carlo Hou", Koudansha, (1994).