分子シミュレーション研究会会誌 “アンサンブル” Vol. 10, No. 4, October 2008 (通巻 44 号)

29

分子動力学シミュレーションにおける温度・圧力制御
第 1 回：能勢の熱浴と能勢・フーバー熱浴
BioMaPS Institute for Quantitative Biology, Rutgers University奥村久士 okumura@biomaps.rutgers.edu
1

はじめに

3 拡張系の方法

今月号から分子動力学シミュレーションにおける温分子動力学シミュレーションにおいて温度一定のア度制御と圧力制御の方法について連載記事を執筆する

ンサンブルや圧力一定のアンサンブルなどミクロカノことになりました．分子動力学シミュレーションにおニカルアンサンブル以外のアンサンブルを生成する場

ける温度制御といえば故能勢修一慶應義塾大学教授の合に (1) 式のハミルトニアン（これを物理系とよぶ）開発された能勢の熱浴が有名です．もし能勢先生がごに熱浴や圧力浴の自由度を追加し温度や圧力を制御す健在でいらっしゃったならば私のような若輩がこの連載る方法がある．物理系と熱浴などの外系をあわせて拡記事を執筆することはなかったと思います．能勢先生

張系とよび，これらの方法を総称して拡張系の方法とが早く亡くなられたことは大変に悲しく残念です．数いう．

十年にわたる研究の歴史があるテーマなので恐縮ではありますが，能勢研究室の先輩で「アンサンブル」編「拡張系」=「物理系」+「外系（例：熱浴）」(5)

集委員でもいらっしゃる森下徹也博士（産総研）に勧拡張系の方法のうち温度を制御してカノニカルアンサめていただきましたので，学生時代に能勢先生から教ンブルを得る代表的な方法が能勢の熱浴である [1, 2]．

えていただいたことをもとに温度・圧力制御の方法について解説したいと思います．私の知識がいたらないところもあるかと思いますが皆様のご指導ご鞭撻のほどよろしく願いいたします．

2

4 能勢の熱浴と能勢・フーバー熱浴
能勢の熱浴では運動量 pi と時間ステップ幅 ∆t を

pi = p′i /s , ∆t = ∆t′ /s(6)

ミクロカノニカルアンサンブルでの分子動力学シのようにスケールし，p′ を仮想運動量，t′ を仮想時間

ミュレーション
分子動力学シミュレーションで最も簡単に得られる
統計アンサンブルはミクロカノニカルアンサンブルである．多くの場合ポテンシャルエネルギー U は座標のみに依存し，そのハミルトニアンは

H0 =N
∑p2i

2mii=1+ U (r)という．そして物理系と熱浴を結びつけた系のハミル
トニアン HN を以下のように導入する．
N
∑p′2
Ps2i
HN =+U(r)++ gkB T0 log s2mi s22Q
i=1(7)(1)

ここで右辺第 1 項と第 2 項は物理系のハミルトニアンとあらわされる．ここで i 番目の粒子の座標，運動量，(1) 式に対応し，第 3 項と第 4 項が熱浴をあらわす．Psは能勢の熱浴の自由度 s に対する正準共役な運動量で質量をそれぞれ r i , pi , mi とあらわし，N 個の粒子の座ある．Q は熱浴の質量に相当し s の時間発展を制御す標の集合を r ≡ {r 1 , r 2 , · · · , r N } と書いた．このハミる．kB はボルツマン定数，T0 は設定温度である．g はルトニアン H0 による運動方程式は正準方程式により系の自由度をあらわすパラメーターで N 個の原子か∂H0pdr i== idt∂pimi∂H0dpi=−= Fidt∂r iらなる系では，仮想時間 t′ について時間発展する場合(2)(3)とあらわされる．F i は i 番目の粒子に働く力で∂UFi ≡ −∂r ig = 3N とする．r, p′ , s, Ps を正準変数として HN から正準方程式を導くと仮想時間 t′ での運動方程式が以下のように決まる．

(4)と定義される．この運動方程式を数値的に解くことによりミクロカノニカルアンサンブルが得られる．

g = 3N + 1，現実の時間 t について時間発展する場合∂HNp′idr i==′′dt∂pi
mi s2dp′i∂HN=−= Fidt′
∂r i

(8)
(9)

30

The Molecular Simulation Society of Japan

ds=dt′

∂HNPs
(10)
=
∂Ps
Q)
(dPs
∂HN1 ∑ p′2i=−− gkB T0=(11)dt′∂ssms2ii
仮想時間 t′ についての運動方程式から現実の時間 t についての運動方程式を求めると

ṙ i =pimiṡṗi = F i − pisPsṡ = sQN∑p2iP˙s =− gkB T0
mii=15 カノニカルアンサンブルが実現することの証明
ハミルトニアン HN はある一定値 E に保たれるので，拡張系全体ではミクロカノニカルアンサンブルになる．そのため拡張系全体の分配関数は∫ ∞ ∫ ∞∫∫
Z =dsdPs dr dp′−∞

0(12)( ())
p′P2× δ H0 r,+ s + gkB T0 log s − E (22)s2Q(13)と書ける．ここで仮想運動量 p′i から現実の運動量 pi(14)(15)への変数変換 pi = p′i /s をおこなう．各粒子の運動量の x, y, z 座標がdpix = dp′ix /s , dpiy = dp′iy /s , dpiz = dp′iz /s (23)のようにスケールされているのでとなる．ここで上付きの点 ( ˙ ) は現実の時間 t での微dp′ = s3N dp

分 d/dt を意味する．(8)−(11) 式および (12)−(15) 式が能勢の熱浴の運動方程式である．

Hoover は現実の時間 t での運動方程式を変形してζ≡

1 dsdsPs= ′ =s dtdtQ(16)とおいた [3]．ζ を用いると能勢の熱浴の運動方程式はpi
miṗi = F i − ζpi(N)1 ∑ p2iζ̇ =− gkB T0
Q i=1 miṙ i =(17)と書ける．この変数変換をおこなうと∫ ∞ ∫ ∞∫∫
Z =dsdPs dr dp s3N0−∞()P2× δ H0 (r, p) + s + gkB T0 log s − E (25)2Qとなる．ここで δ 関数の恒等式

δ(s − s0 )|f ′ (s)|(18)δ(f (s)) =(19)を用いる．今の場合{(の 3 式で閉じた形になる．これを能勢・フーバー熱浴(24)

s0 = exp

ただし

P2E − H0 (r, p) − s2Q

f (s0 ) = 0(26)})/gkB T0(27)という．
瞬間温度はT (t) =N1 ∑ p2i
gkB i=1 mi(20)とあらわせるのでこれを使って (19) 式を書きかえると

dζ
gkB
=
(T (t) − T0 )dt
Q(21)となる．瞬間温度 T (t) が設定温度 T0 より高ければ ζは大きくなる．(18) 式で ζ は抵抗係数のような役割を果たしているので，粒子の運動量 pi は遅くなる．ただであるから δ 関数の積分を実行すると∫ ∞∫∫1Z =dPs dr dp
gkB T0−∞{()}
3N + 1P2× exp
E − H0 (r, p) − s(28)
gkB T02Qとなる．ここで g = 3N + 1 とすると()
∫ ∞1E − Ps2 /2Q
Z =dPsexp
gkB T0kB T0−∞}{∫∫
H0 (r, p)
× dr dp exp −kB T0{z}|(29)し抵抗係数とはいってもこの ζ は正の値も負の値も両カノニカルアンサンブル方とりうることに注意する．逆に瞬間温度 T (t) が設定であるから仮想時間 t′ でサンプルする場合には g =温度 T0 より低くければ ζ は小さくなり，粒子の運動づくように運動量にフィードバックがかかっているこ

3N + 1 とすれば座標と運動量の両方にについてカノニカルアンサンブルが得られることがわかる．
一方，現実の時間 t でサンプルする場合には g = 3Nとがわかる．とするとカノニカルアンサンブルが得られる．このこ量 pi を速くする．瞬間温度 T (t) が設定温度 T0 に近分子シミュレーション研究会会誌 “アンサンブル” Vol. 10, No. 4, October 2008 (通巻 44 号)とを次に示す．現実の時間 t でサンプルする場合，物理量 A(r, p) の平均値は∫
1 τlimdt A (r(t), p(t))τ →∞ τ 0∫ τ1= lim
dt A (r(t), p′ (t)/s(t))τ →∞ τ 0∫ ′τ ′ 1 τ ′ A(r(t′ ), p′ (t′ )/s(t′ ))= limdt
τ →∞ τ τ ′ 0s(t′ )実際に能勢・フーバー熱浴をもちいて分子動力学シ数値積分しなければならない．能勢・フーバー熱浴の運動方程式を差分化する方法にはいくつかあるが，そのうちプログラムを書くのがもっとも簡単なのは Martyna
らが提案した時間反転可逆な積分法であろう [4]．この方法を次に説明する．
能勢・フーバー熱浴における時間発展は (17)−(19) 式と計算される．ここで dt = dt /s であるので τ は

∫ τ
τ=

∫ τ
dt =

001 ′dtsと表される．これを用いると
∫
1 τ
lim
dt A (r(t), p(t))τ →∞ τ 0∫ ′1 τ ′ A(r(t′ ), p′ (t′ )/s(t′ ))limdtτ ′ →∞ τ ′ 0s(t′ )
=
∫ τ′
11
lim
dt′
τ ′ →∞ τ ′ 0s(t′ )⟨A(r, p′ /s)/s⟩=⟨1/s⟩

6 能勢・フーバー熱浴の時間発展
ミュレーションをおこなうには運動方程式を差分化し′′31で記述されるので位相空間は (r, p, ζ) で張られる．よっ(30)て物理量 A(r, p, ζ) の時間発展はȦ(r, p, ζ) =∑ir˙i ·∂A ∑∂A ∑ ∂Aṗi ·ζ̇++(34)∂r i∂pi∂ζiiと書ける．ここで演算子

D≡∑ir˙i ·∑∑ ∂∂∂+ṗi ·+ζ̇∂r i∂pi∂ζii(35)

を導入すると (34) 式は(31)となる．ここで，⟨· · ·⟩ は仮想時間 t′ での平均を意味すȦ(r, p, ζ) = DA(36)と書き換えられる．この微分方程式の形式解はる．⟨A(r, p′ /s)/s⟩ および ⟨1/s⟩ を計算する際 (25) 式

A(t + ∆t) = eD∆t A(t)の s3N が s3N −1 になることに注意すると，(28) 式の計(37)となる．ここで eD∆t を時間発展演算子とよぶ．この算と同様に⟨A(r, p′ /s)/s⟩
⟨1/s⟩
{}
∫
∫
3N
dr dp A(r, p) exp −H0 (r, p)
gkB T0{}∫∫3N + 1
dr dp exp −H0 (r, p)
gkB T0{}∫∫=3N
dr dp exp −H0 (r, p){ gkB T0}∫∫
3N + 1
dr dp exp −H0 (r, p)
gkB T0{}∫∫3N
dr dp A(r, p) exp −H0 (r, p)
gkB T0{}∫∫=3N
dr dp exp −H0 (r, p)
gkB T0(32)

形式解は厳密に正しいがこのままでは数値積分ができない．そこで (17)−(19) 式を (35) 式に代入し以下のように演算子 D を分割する．

D = D1 + D2 + D3(38)
(N)
2∑ p
∑
∂
1pi
∂
i
D1 =·+− gkB T0mi ∂r i
Q i=1 mi∂ζi(39)∑∂D2 =Fi ·∂pii∑∂
D3 = −ζpi ·∂pii

(40)
(41)

演算子 D を分割したのにともない，時間発展演算子(33)eD∆t も近似的に以下のように分割する．これを鈴木・トロッター分割という．
()∆t∆t∆t∆t
3eD∆t = eD3 2 eD2 2 eD1 ∆t eD2 2 eD3 2 + O (∆t)となる．ここで，⟨· · ·⟩N V T はカノニカルアンサンブル(42)
この分割の仕方は他にも考えられるがここでは最も簡での平均である．すなわち，現実の時間 t でサンプル単な分割法を用いる．eD1 ∆t はを得られる．ゆえに g = 3N とすれば
∫1 τdt A (r(t), p(t)) = ⟨A(r, p)⟩N V Tlim
τ →∞ τ 0

する時は g = 3N とすればカノニカルアンサンブルが得られる．

eD1 ∆t = 1 + D1 ∆t +1 2D (∆t)2 + · · ·2! 1(43)32

The Molecular Simulation Society of Japan
2.0と展開されるので r i , ζ に演算するとiQi=1mi1.50.8)− gkB T01.0Ts(44)0.60.5Q = 0.10.4∆t (45)Q = 0.10.0
1.01.5となる．ここで D1 は p に演算すると 0 なので ∆t の0.81.0TeD1 ∆t ζ = ζ +pi∆tmi(N1 ∑ p2seD1 ∆t r i = r i +1.00.62 乗以上の高次項は 0 になることに注意する．つまり(44), (45) 式は厳密に正しい．同様に

0.5

Q=10.4

Q=10.01.01.5(46)0.8となる．D3 の演算は ∆t の高次項が 0 にならないが，0.6

公式

0.40.5Q = 10Q = 100.01.0(47)1.50.8sを用いると1 2D (∆t)2 + · · ·}pi
2! 3
1= pi {1 + (−ζ∆t) + (−ζ∆t)2 + · · ·}2!
= pi e−ζ∆t(48)

eD3 ∆t pi = {1 + D3 ∆t +1.0T1
e = 1 + x + x2 + · · ·2!
x1.0TseD2 ∆t pi = pi + F i ∆t

0.60.5Q = 100

0.4

Q = 100

0.00

123time450

123time45

図 1: 能勢の熱浴の自由度 s と瞬間温度 T の時間変化．と計算される．各時間発展演算子による時間発展が計算できたので (42) 式の順番で r, p, ζ を時間発展させ粒子間相互作用 u(rij ) にはレナード・ジョーンズポテ

る．実際に分子動力学シミュレーションをおこなうにンシャルを用いた．

は以下の手順でシミュレーションプログラムを書けばよい．

∆tpi ← pi + F i(2∑u(rij )i>j()∆tpi ← pi exp −ζ2∆tpi ← pi + F i
2pi
ri ← ri +∆tmi(N)1 ∑ p2i− gkB T0 ∆tζ ← ζ+Q i=1 mi

U (r) ={((49)u(rij ) = 4ϵσrij(55))12(−σrij)6 }
(56)

(50)

ここで rij = |r i − r j | は 2 粒子間距離，ϵ はポテンシャ(51)

ルエネルギーの深さ，σ はレナード・ジョーンズ粒子の直径である．今後は ϵ と σ を用いてそれぞれエネル(52)ギーと長さを無次元化した単位系を用いる．周期境界条件を用い，相互作用のカットオフ距離は rc = 4.0，粒(53)子数は N = 500，密度は ρ = 0.8 に設定した．設定温(54)の密度と設定温度ではレナード・ジョーンズ粒子系は

ここで左向きの矢印 ← はプログラム中での代入を意味液体になる [5]．熱浴の質量は Q = 0.1 から Q = 100pi ← pi exp −ζ∆t2)

度は T0 = 1.0，初期温度は T (t = 0) = 0.5 とした．こする．(49)−(54) 式のプロセスを 1 回計算することにより能勢・フーバー熱浴を用いた分子動力学シミュレーションを 1 ステップおこなうことができる．(49)−(54)式において熱浴がない，すなわち ζ = 0 の場合には速度ベルレ法によるミクロカノニカルアンサンブルでの時間発展の表式が得られる．

まで 4 つの値に設定した．時間ステップ幅 ∆t = 0.002としてシミュレーションを 2500 ステップおこなった．
能勢の熱浴の自由度 s と瞬間温度 T の時間変化を図1 に示す．Q が小さくなるほど s の振動は速く，Q が大きくなるほど s の振動は遅くなる．Q = 0.1 や Q = 100の場合，t = 5 までの範囲では温度が平衡に達していない．一般に s の振動は速すぎても遅すぎても熱浴と応用例：レナード・ジョーンズ流体

物理系の間でエネルギーのやり取りが効率的におきず，前節で述べた時間発展法を用いて実際に分子動力学平衡に達するまでに時間がかかってしまう．s の振動のシミュレーションをおこなった例について述べる．2 体

時間スケールが粒子の運動の時間スケールと同程度で7分子シミュレーション研究会会誌 “アンサンブル” Vol. 10, No. 4, October 2008 (通巻 44 号)

あるときに速く平衡に達する．この例では Q = 1 ∼ 10程度に設定するのがよい．Q = 1 と Q = 10 の場合ではそれぞれ t = 0.5 と t = 1.0 程度で平衡に達し瞬間温度 T は設定温度 T0 = 1.0 の周りで小刻みに揺らいでいる．

8

まとめ
この解説記事では能勢の熱浴と能勢・フーバー熱浴について説明した．その運動方程式を導出し，座標と運動量の両方についてカノニカルアンサンブルが得られることを証明した．さらに能勢・フーバー熱浴を用いた分子動力学シミュレーションを実行するために運動方程式を差分化し数値積分する手順も紹介した．またレナード・ジョーンズ流体の分子動力学シミュレーションをおこない熱浴の質量 Q の値の選び方を述べた．この解説記事を読んでいただければ能勢の熱浴や能勢・フーバー熱浴の運動方程式を理解できるだけでなく，初学者でも自分でシミュレーションプログラムを書いて実行できるようになると思う．次回以降も分子動力学シミュレーションの経験があまりない人でもプログラムが書けるように，なるべくわかりやすく解説していきたい．

参考文献[1] S. Nosé, Mol. Phys., 52, 255 (1984).
[2] S. Nosé, J. Chem. Phys., 81, 511 (1984).
[3] W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
[4] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M.
L. Klein, Mol. Phys., 87, 1117 (1996).
[5] H. Okumura and F. Yonezawa, J. Chem. Phys., 113,9162 (2000).
著者紹介

奥村久士 (博士 (理学))：〔経歴〕

2002 年慶應義塾大学理工学研究科博士課程修了，東京大学，分子科学研究所，名古屋大学を経て 2008 年から現所属．
〔専門〕統計力学，分子シミュレーション．
〔趣味〕海外旅行，落語鑑賞．

33