分子シミュレーション研究会会誌 “アンサンブル” Vol. 11, No. 2, April 2009 (通巻 46 号)

43

分子動力学シミュレーションにおける温度・圧力制御
第 3 回：速度スケーリング法，ガウス束縛法，ベレンゼン熱浴
自然科学研究機構
分子科学研究所奥村久士

1

はじめに前回，前々回は温度を制御する拡張系の方法として能勢の熱浴，能勢・フーバー熱浴，能勢・ポアンカレ
熱浴の解説をした．今回は拡張系の方法以外で温度を制御する方法である速度スケーリング法，ガウス束縛法，ベレンゼン熱浴について解説する．

2)(∆t∆tr i (t + ∆t) = r i (t) + v i t +2

レーションを行う場合でも，その前の準備段階でこの分子動力学シミュレーションで温度制御をおこなう

3 ガウス束縛法

もっとも簡単な方法は速度スケーリング法である [1]．
速度スケーリング法では各ステップごとに運動エネル
ギーが指定した温度に対応するように速度をスケール

3.1 ガウス束縛法の運動方程式
温度を制御するため，運動エネルギー K にする．例としてリープ・フロッグ法で時間発展させる

K=場合について述べる．リープ・フロッグ法では，粒子i に働く力 F i と粒子の質量 mi を用いてまず速度を()()∆t∆t
F i (t)′vi t +(1)= vi t −+∆t
22mi()と進ませる．速度 v ′i t + ∆tに対応する運動エネル
2′iN
∑
p2i2m
i=1

= 一定(8)の束縛を加えて運動方程式を解く方法もある．これは温度一定の ガウス束縛法と呼ばれている [2, 3, 4]．
今後，3N 次元のベクトルを考え下付き添字のない運動量 p，座標 r ，粒子間力 F をそれぞれ，p =(p1 , p2 , · · · , p3N )T , r = (r1 , r2 , · · · , r3N )T , F =ギー K は()∆t3Nv ′2t+=kB T ′i222(7)と座標を時間発展させる．他の温度一定の方法でシミュ方法を用いることが多い．

∑ mihokumura@ims.ac.jp

として (3) 式のスケーリングを行えばよい．この速度()v i t + ∆tを用いて2速度スケーリング法

K′ =計算科学研究センター

理論・計算分子科学研究領域(2)(F1 , F2 , · · · , F3N )T とする．ここで上付き添字の T は転置を意味する．すると，運動エネルギー K はである．ここで，T ′ は運動エネルギー K ′ に対応する

K=′瞬間温度である．この T が指定する温度 T0 になるよ

p2= 一定2m(9)うに，すなわち運動エネルギー K ′ が K = 3N
2 kB T0 にと書ける．ここで束縛力が最小になるように束縛力をなるように速度を s 倍スケールする：
()
()
∆t
∆t
vi t += sv ′i t +.
22かける．一般に束縛条件は(3)

R(r, ṙ, t) = 0

s の値は瞬間温度 T ′ から求まる．運動エネルギーは)()∑ mi (∑ mi∆t∆t22′2(4)v t+=svt+2 i22 i2iiと書ける．この時間微分をとると

dR
∂R ∂R
∂R=+· ṙ +· r̈ = 0 .
dt
∂t∂r
∂ ṙ

(10)

(11)

束縛力 F c がかかった場合の運動方程式

であるから

3N
3NkB T0 = s2
kB T ′ .
22
よって√s=T0T′mr̈ = F + F c(12)(5)の両辺に ∂R/∂ ṙ をかけて (11) 式を代入すると

(6)

Fc ·

∂R
∂R∂R
∂R= −F ·−m− mṙ ·∂ ṙ
∂ ṙ
∂t∂r

(13)

44

The Molecular Simulation Society of Japan

よって未定係数 ζ が.
rN∑p

束縛力 Fc

原子間力 Fi

ζ=−mi=1N∑
p· Fi=−N∑
p2· Fi(22)

gkB T0ii=1imi=1mと求まる．(17), (18), (22) 式がガウス束縛法の運動方程式である．
図 1: ガウス束縛法の模式図.

力 F i はポテンシャルエネルギー U の座標微分となる．ここで右辺は粒子の座標と運動量 (r, ṙ) が決Fi = −まれば決まる値であり，束縛力 F c の選び方にはよらであるので，(22) 式の分子はない．束縛力の絶対値 |F c | を最小にするには∂RFc ∝∂ ṙN∑pmi=1(14)とすればよい．すなわち束縛条件 R の面に垂直にとればよい（図 1 参照）．

m 2 gṙ − kB T0 = 022アンサンブルが得られることがわかる．(15) 式の ṙ 微分を計算すると位相空間 Γ = (r, p) での状態分布関数 f (Γ ) を考える．状態分布関数 f (Γ ) について連続の式∂f∂+· (Γ̇ f ) = 0∂t∂Γg− kB T0 = 02う．(26) 式左辺第 2 項は

imi=1· ṗi = 0 .

imi=1· (F i − ζpi ) = 0 .

(∂· Γ̇∂Γ)
f(27)と書きかえれらる．(26), (27) 式を用いると f (Γ ) の全(18)

微分はdf
∂f
∂f=+ Γ̇ ·=−dt
∂t∂Γ
(19)

(20)

これに (18) 式を代入すると
N∑p∂∂f· (Γ̇ f ) = Γ̇ ·+∂Γ
∂Γ

(17)

を満たすように決める．(19) 式の時間微分をとって
N
∑
p(26)

が成り立つ．これを一般化されたリウビル方程式といと書ける．未定係数 ζ は温度一定の束縛条件
i(25)3.2 ガウス束縛法により得られる統計アンサンブルばよい．そこで ζ を未定係数として運動方程式は2mi=1(24)dU
ζ = − dt .
gkB T0(16)なので，束縛力 F c を運動量 p に比例するようにとれN∑p2i=1とも書ける．

論から g = 3N − 1 とすれば座標空間でのカノニカルdr i
p= idtmdpi= F i − ζpidt∂UdU=∂r idtr˙i ·時間全微分を用いて(15)と書ける．ここで g は系の自由度をあらわし，後の議∂R= mṙ = p∂ ṙ

· Fi =N∑(23)と書ける．すなわち ζ は全ポテンシャルエネルギーの温度一定の束縛がかかっている場合，束縛条件 R はR(r, ṙ, t) =i∂U∂r i(21)(∂· Γ̇∂Γ)f(28)と書ける．(28) 式右辺の ( ) 内は (17), (18) 式を用いて以下のように計算される．
)N (∑∂∂∂· Γ̇ =· ṗ +· ṙ i∂Γ∂pi i ∂r ii=1N {∑∂∂=· Fi −· (ζpi )∂p∂piii=1}∂ ( pi )+·∂r im= −3N ζ −N∑

pi ·i=1= −(3N − 1)ζ .

∂ζ∂pi(29)分子シミュレーション研究会会誌 “アンサンブル” Vol. 11, No. 2, April 2009 (通巻 46 号)

45

ここで 2 行目の計算には(37)ṗi = F i∂
· F i = 0,
∂pi

∂
·p =0
∂r i i(30)の代わりにベレンゼン熱浴ではを，3 行目の計算には

 N
∑ pj· Fj NNm∑∑∂ζ∂ 
 j=1 jpi ·= −pi · N2∂p∂p∑
iipj i=1i=1
m
jj=1N∑piṙ i =ṗi = F i +(39)

T が設定温度 T0 にどの程度速く近づくかを制御する．

K=∑ p2i2mii

mj
j=1()
NN
∑∑
pj
p2i

· Fjmjmij=1i=1+2
2N
∑p2j

mj
j=1=3NkB T2(40)の関係がある．(39) 式右辺第 2 項の由来は以下のように解釈できる．瞬間温度 T の時間変化はṪ =(31)を用いた．(25), (28), (29) 式より(32)であるから g = 3N − 1 と選べば

ゆえに

12τ)
T0− 1 piTという方程式を使う．ここで τ は緩和時間で瞬間温度

mi= − i=1N
∑ p2j

3N − 1 1 dUdf=−fdtg
kB T0 dt

(38)
(T は運動エネルギー K と· Fi= ζ − 2ζ = −ζpimi

df1 dU=−f.
dt
kB T0 dt

(33){}
U (r)f (r, p) = C exp −.
kB T0(34)2 ∑ pi · ṗi
3N kB imi()∑ 2pi
2 ∑ pi2
1 T0=· Fi +−1
3N kB i mi
3N kB 2τ Tmii
2 ∑ pi1=· F i − (T − T0 )(41)
3N kB i miτと計算される．(41) 式右辺第 1 項はミクロカノニカルMD でもあらわれる項である．(41) 式右辺第 2 項がベレンゼン熱浴の特徴で，瞬間温度 T が設定温度 T0 に近づくように( )dT
1= − (T − T0 )
(42)
dt 熱浴
τという熱浴をつけたことに相当している．

ここで C は定数である．最後に，運動エネルギー一定(38), (39) 式を時間発展させるには予測子・修正子法の束縛があることから状態分布関数 f (r, p) は(N){}
∑ p2g
U (r)
if (r, p) = δ− kB T0 exp −2mi2kB T0i=1(35)などの方法を用いてもよいが，以下のように近似的にと決まる．ガウス束縛法では座標空間 r ではカノニカルアンサンブルが得られるが，運動量空間 p はカノニカルアンサンブルにならない．しかし，実際の MD シミュレーションでは座標空間でカノニカルアンサンブ

速度スケーリング法を用いることが多い．例えば，リープ・フロッグ法を基に時間発展させる場合にはまず運動量を(())∆t∆tt+= pi t −+ F i (t)∆t(43)22)(と進ませる．運動量 p′i t + ∆tによる瞬間温度2p′i

T =ルが得られていれば問題ないことが多く，ガウス束縛法も能勢・フーバー系の方法などと同じように良く使われている．

4ベレンゼン熱浴
最後にベレンゼン熱浴について説明する．ミクロカ

pṙ i = imiを求め，スケーリングファクター s を()∆t T0−1s=1+2τ T
またはノニカルアンサンブルにおける運動方程式(36)2 ∑ p2i
3N kB i 2mi()] 21[∆t T0−1s= 1+
τT

(44)

(45)

(46)

46

The Molecular Simulation Society of Japanと計算する．このスケーリングファクター s を用いて運動量を s 倍スケールする．
()()∆t∆t′pi t += spi t +.
2
2[1] L. V. Woodcock, Chem. Phys. Lett. 10 (1975) 257.
(47)

最後にスケールした運動量 pi を用いて座標を時間発展させる．

(r i (t + ∆t) = r i (t) +[2] W. G. Hoover, A. J. C. Ladd, and B. Moran, Phys.
Rev. Lett. 48 (1982) 1818.
[3] D. J. Evans, J. Chem. Phys. 78 (1983) 3297.

)∆t

pi t + 2
∆t .
mi

参考文献

(48)

[4] S. Nosé, Prog. Theor. Phys. Suppl., 103 (1991) 1.
[5] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gun-(43)−(48) 式の時間発展が (38), (39) 式の運動方程式に対応していることは以下の議論からわかる．(47) 式に (43) 式と (45) 式を代入すると
()
{()}
∆t
∆t T0pi t += 1+−122τ T{ ()}∆t× pi t −+ F i (t)∆t2()∆t= pi t −+ F i (t)∆t2() ()1 T0∆t+− 1 pi t −∆t2τ T2

steren, A. DiNola, and J. R. Haak, J. Chem. Phys., 81,3684 (1984).
[6] T. Morishita, J. Chem. Phys., 113, 2976 (2000).
奥村久士 (博士 (理学))：〔経歴〕

2002 年慶應義塾大学理工学研究
著者紹介科博士課程修了，東京大学日本学術振興会特別研究員 (PD)，分子科学研究所助手，名古屋大学(49)

COE 特任講師，ラトガース大学研究助手，同研究助教授を経てとなる．これは ∆t のオーダーまで (39) 式と一致して現在，自然科学研究機構計算科いる．

学研究センター准教授〔専門〕統+O((∆t)2 )ベレンゼン熱浴により生成されるアンサンブルはカノニカルアンサンブルでもミクロカノニカルアンサンブルでもない．そのアンサンブルは緩和時間 τ に依存する．座標空間について τ → 0 の極限でカノニカルアンサンブルに，τ → ∞ の極限でミクロカノニカルアンサンブルに近づくことが森下により示されている [6]．

5まとめ
今回は拡張系の方法以外で温度を制御する分子動力学法（速度スケーリング法，ガウス束縛法，ベレンゼン熱浴）を解説した．一般にシミュレーションの初期段階では速度スケーリング法を用いて温度制御がおこなわれる．系が平衡に達した後，データをサンプルするには座標空間，運動量空間ともにカノニカルアンサンブルを生成できる能勢の熱浴，能勢・フーバー熱浴，または能勢・ポアンカレ熱浴を用いるのがよい．したがって，これから自分でプログラムを書こうという人はガウス束縛法やベレンゼン熱浴を勉強しても使わないかもしれない．しかし，これらの手法はともに分子シミュレーションの発展の歴史において重要な役割を果たしてきた．また市販または配布されているシミュ
レーションプログラムの中には能勢の熱浴などが入っていないものもある．その場合にはガウス束縛法かベレンゼン熱浴を使うことになると思うので，必要に応じて今回の解説を読んでいただければ幸いである．

計力学，分子シミュレーション．
〔趣味〕海外旅行，落語鑑賞．