Journal of Computer Chemistry, Japan
Online ISSN : 1347-3824
Print ISSN : 1347-1767
ISSN-L : 1347-1767
技術論文
Mathematicaによる分子動力学
山田 祐理片岡 洋右
著者情報
ジャーナル フリー HTML
電子付録

2015 年 14 巻 5 号 p. 172-176

詳細
Abstract

分子動力学(MD)法の初学者への入門のために,Wolfram MathematicaによるMDのプログラムを開発した.MD計算やその結果を解析するアルゴリズムを理解し易いように,扱う分子間相互作用はLennard-Jones 12–6ポテンシャル,MDセルは立方体,アンサンブルはNEVまたはNVTアンサンブルのみとした.Mathematicaのノートブックファイル内で,粒子数,数密度,温度などの計算条件を設定する.計算結果は,熱力学量のほかに粒子の軌跡,粒子配置,二体相関関数,速度自己相関関数,平均自乗変位および自己拡散係数が出力される.

1 はじめに

Wolfram Mathematica [1]は元来「数式処理」ソフトウェアであるが,非常に多くの機能を備えており,数学のみならず自然科学全般において幅広く活用することができる.また,豊富な組み込み関数を持ち,プログラミング言語としても洗練されている.

我々はこのMathematicaを用い,分子動力学(MD)法をこれから学ぶ初学者を想定して,MDの入門プログラムを開発した.これは,以前に報告したExcel VBAによるMDプログラム [2,3]と,目的も内容も共通している.Mathematicaによるプログラムは,入出力も含めすべて一つのファイルにまとめられ,また組み込み関数を活用することで,計算のアルゴリズムも分かり易く示すことができる.本プログラムが,MDプログラムの基本的構造や,計算および解析の手法を学ぶ一助になることを期待する.

2 モデル

本プログラムでは簡便のために,Lennard-Jones (LJ) 12–6ポテンシャル   

u i j ( r i j ) = 4 ε { ( σ r i j ) 12 ( σ r i j ) 6 } (1)
に従う一種類の粒子からなる系のみを扱う.riji番目とj番目の粒子間の距離,εσはそれぞれエネルギーと長さの次元を持つパラメータである.本プログラムにおいて,各種物理量はこのεσ,ボルツマン定数k,およびLJ粒子の質量mによって無次元化された値を用いる:    エネルギー/ε,長さ/σ,温度/εk−1       圧力/εσ−3,数密度/σ−3,時間/τ    (2)    ただしτ = σ(m/ε)1/2である.

3 MDプログラム

プログラムの開発にはMathematica 9.0.1.0を用いた.プログラムの基本構造は成書(片岡,1994) [4]を元にしており,Mathematicaの仕様に合わせて修正を加えた.

3.1 初期配置と計算条件

MDセルの形は立方体に限定する.まず,指定した数密度numberDensity※に対応した面心立方(fcc)構造の単位格子を用意し,これをxyzの各方向にそれぞれnfcc個積み重ねる.したがって,MDセルに含まれる粒子数nspは,次のように定まる.    nsp = 4 × nfcc3.    (3)   

こうして作ったfcc構造に対して,最近接粒子間距離(21/6 × numberDensity−1/3)のdelta倍を上限幅として,全粒子をそれぞれランダムに変位させた構造を,MDの初期配置とする.Figure 1に初期配置の例を示した.

Figure 1.

 An example of the initial configuration at numberDensity = 0.6, nfcc = 4 (number of particles nsp = 4 × nfcc3 = 256), and delta = 0.1.

計算条件の設定は,ノートブック内の設定値を書き換えることで行われる.設定できる量は,numberDensitynfccdeltaのほか,温度(NEV計算の場合は初期温度) temp,温度制御をするかどうかのフラグtconst,MDステップ数nstep,MD時間刻みdtなどである.

※  このフォントで表された量は,Mathematicaノートブック内で使用する変数であることを示す.

3.2 計算方法

いずれもMDのごく基礎的な方法であるので,詳しくは成書[5, 6など]もしくは既報 [7]を参照されたい.

運動方程式の数値積分にはVerlet法を用い,NVT計算における温度制御には速度スケーリング法を用いた.本プログラムでは,NVT計算の全ステップで速度スケーリングを行っている.MDセルには三次元周期境界条件およびminimum image conventionを適用する.相互作用のカットオフ距離は,立方体セルの幅の半分とした.

4 プログラムの使用法

Mathematicaがインストールされたコンピュータであれば,付録のノートブックファイル [8]を開いて「評価 > ノートブックを評価」すれば,既定の計算条件numberDensity = 0.6,nfcc = 4 (nsp = 256), temp = 1,nstep = 1000,dt = 0.01,delta = 0.1におけるNEV計算が実行できる.計算の途中経過および結果は,ノートブック内「5. 分子動力学シミュレーション」および「6. 分子動力学シミュレーション結果の解析」に順次出力される.計算条件を変えるには,ノートブック内「3. 計算条件の設定」にて各設定値を書き換えればよい.

Mathematicaのない環境でも,無料のWolfram CDF Player [9]をインストールすれば,閲覧専用ではあるがノートブックファイルを開くことができる.

5 液体の計算例

液体とみなせる状態の計算例を以下に示す.設定値numberDensity = 0.6,nfcc = 4 (nsp = 256), temp = 1,nstep = 1000,dt = 0.01,delta = 0.1は共通とし,NEVNVTの計算をそれぞれ行った.この条件での計算にかかる実時間は,現在の一般的な家庭用PCで15∼30分程度である.

5.1 熱力学量の時間発展;NEVNVTの比較

NEV計算とNVT計算の比較のため,一粒子あたりのハミルトニアンおよび系の温度の時間発展をFigure 2に示す.NEV計算ではハミルトニアンが,NVT計算では温度が保存されている.プログラム中では,これらの他に運動エネルギー,ポテンシャルエネルギー,圧力などが出力される.

Figure 2.

 Time development of Hamiltonian per particle and temperature of the system at numberDensity = 0.6, nfcc = 4, temp = 1, nstep = 1000, dt = 0.01, and delta = 0.1.

5.2 粒子の軌跡

以下,結果はすべてNEV計算のものを示す.上記の設定値においては,NVT計算の結果も大きくは変わらない.

粒子のx座標の変化の例をFigure 3に,粒子の軌跡の例をFigure 4にそれぞれ示した.周期境界条件により,セルからはみ出した粒子がセル内に戻されているのが分かる.

Figure 3.

 Time development of x-coordinate of a particle at numberDensity = 0.6, nfcc = 4, temp = 1, nstep = 1000, dt = 0.01, and delta = 0.1.

Figure 4.

 Trajectories of some particles at numberDensity = 0.6, nfcc = 4, temp = 1, nstep = 1000, dt = 0.01, and delta = 0.1.

5.3 構造

Figure 5に最終配置の例を示す.初期配置(Figure 1)に見られるような規則的構造は見て取れず,液体的であることが分かる.

Figure 5.

 An example of the final configuration at numberDensity = 0.6, nfcc = 4, temp = 1, nstep = 1000, dt = 0.01, and delta = 0.1.

Figure 6は二体相関関数である.比較のために初期配置の二体相関関数(点線)も示した.この図からも,初期配置における規則的構造が崩れていることが分かる.

Figure 6.

 Pair correlation function at numberDensity = 0.6, nfcc = 4, temp = 1, nstep = 1000, dt = 0.01, and delta = 0.1. Dashed curve indicates the pair correlation function of the initial configuration.

5.4 自己拡散係数

MD結果の解析の例として,速度自己相関関数(velocity autocorrelation function; VAF)および平均自乗変位(mean square displacement; MSD)から,それぞれ自己拡散係数Dを求めた.   

D = 1 3 0 v i ( t ) v i ( 0 ) d t = lim t 1 6 t | r i ( t ) r i ( 0 ) | 2 . 4
()

VAFおよびMSDのプロットをFigure 7に示した.(4)に従って両方の方法でDを求めると,次のように互いに良く一致した.    VAFより: D = 0.127σ2τ−1,       MSDより: D = 0.120σ2τ−1. ()    5    この値も,Figure 6と同様に液体の状態を示唆する.

Figure 7.

 Velocity autocorrelation function (VAF) and mean square displacement (MSD) at numberDensity = 0.6, nfcc = 4, temp = 1, nstep = 1000, dt = 0.01, and delta = 0.1.

6 固体および気体の計算について

数密度numberDensityおよび温度(NEV計算では初期温度) tempの設定により,固体および気体に相当する状態点におけるMDも可能である.ただし気体の計算においては,充分に乱雑な初期配置を用意したほうが短時間で緩和するので,deltaは大きめに設定するとよい.Figure 8に,気体に相当する密度における初期配置の例を示した.

Figure 8.

 An example of the initial configuration at numberDensity = 0.001, nfcc = 4, and delta = 0.3.

7 結論

MDの初学者への入門のために,Mathematicaを用いた基礎的なMDのプログラムを開発した.Mathematicaの機能を用いて,相互作用関数の定義,計算条件と初期配置の設定,数値積分の方法,結果の出力と解析の方法などを簡潔なプログラムコードで示した.また,さまざまな状態点におけるMDの結果を比較することによって,固液気の三態における分子の挙動や静的・動的構造の違いが容易く理解できる.

References
  • 1   Mathematica, Wolfram Research, Inc., Champaign, IL (1988–2015).
  • 2   Y. Kataoka, Y. Yamada, J. Comput. Chem. Jpn., 14, 18 (2015).
  • 3   Y. Kataoka, Y. Yamada, Bull. Res. Cent. Comp. Multimedia Stud. Hosei Univ., 29, 1 (2015).
  • 4   Y. Kataoka, "Bunshi Dourikigaku Hou To Monte Carlo Hou," Koudan-sha Scientific (1994).
  • 5   A. Ueda, "Computer Simulation," Asakura Shoten (1990).
  • 6   S. Okazaki, N. Yoshii, "Computer Simulation No Kiso [Dai 2 Han]," Kagaku Dojin (2011).
  • 7   Y. Yamada, Y. Kataoka, Bull. Res. Cent. Comp. Multimedia Stud. Hosei Univ., 29, 22 (2015).
  • 8   3DLJMD.nb
  • 9   http://www.wolfram.com/cdf-player/
 
© 2015 日本コンピュータ化学会
feedback
Top