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

35

分子動力学シミュレーションにおける温度・圧力制御
第 2 回：シンプレクティック分子動力学法と能勢・ポアンカレ熱浴
BioMaPS Institute for Quantitative Biology, Rutgers University
奥村久士 okumura@biomaps.rutgers.edu
1

はじめにの i 成分を Γ i と書くことにする．

今回はシンプレクティック分子動力学法について解説する [1]．シンプレクティック分子動力学法を用いるとシミュレーションを長時間安定におこなうことができる．実時間における能勢の熱浴のシンプレクティッ
ク形式である能勢・ポアンカレ熱浴についても解説すΓ˙i′ =

6N∑∂Γ ′i∂Γjj=1Γ˙j =6N∑

6N∑Mij Γ˙j =j=1

Mij Jjkj,k=16N∑∂Γ ′ ∂H
∂H=Mij Jjk l=
Mij Jjk Mlk∂Γk ∂Γl′∂Γl′
j,k,l=1j,k,l=16N∑(7)

る [2, 3]．

2∂H∂Γkシンプレクティック条件

すなわち′Γ̇ = M J M T

一般座標 qi と一般運動量 pi の組を()qΓ =p∂H∂Γ ′(8)である．ここで M は Γ から Γ ′ への正準変換におけ(1)るヤコビ行列で，その (i, j) 成分はと書くことにする．ここで，q と p はそれぞれ 3N

Mij =個の成分をもっている：q = (q1 , q2 , · · · , q3N ) ，p =T∂Γi′∂Γj(9)(p1 , p2 , · · · , p3N )T ．qi と pi の時間発展は正準方程式により
dqi∂Hdpi∂H=,=−(2)dt∂pidt∂qiで与えられる．(6) 式と (8) 式を比較するとと書ける．ここで H はハミルトニアンである．Γ を用が成り立つ．(10) 式をシンプレクティック条件という．

いてこの時間発展を記述すると

MJMT = J(10)

ある時刻 t における一般座標 q(t), 一般運動量 p(t) がハミルトニアンにしたがって時間発展し，時間 ∆t 後Γ̇ = J∂H∂Γ(3)となる．J は 6N × 6N 行列で()0 1
J=−1 0の時間発展は変数変換
()
q(t)′

Γ =q ′ (q, p)p′ (q, p))ら正準変換により導かれるので Γ ′ についても正準方程式∂HΓ̇ = J∂Γ ′(→q(t + ∆t)p(t + ∆t))(11)とみなすことができる．変数 q(t), p(t) の組も q(t+∆t),p(t + ∆t) の組もともに同じハミルトニアンによる正準方程式にしたがうので，この変換は正準変換である．
上述のように正準変換においてはシンプレクティック(5)への時間に依存しない正準変換を考える．Γ ′ は Γ か′p(t)(4)と定義される．1 は 3N × 3N 単位行列である．Γ から別の座標と運動量の組(にそれぞれ q(t + ∆t), p(t + ∆t) に変化したとする．こ条件が成り立っているので，ハミルトニアンにしたがう時間発展で得られる一般座標，一般運動量はシンプ
レクティック条件 (10) 式を満たしている．現実の物体はシンプレクティック条件を満たしながら運動しているが，運動方程式を差分化し数値積分をおこなうと必(6)ずしもシンプレクティック条件を満たすとは限らない．

が成り立つ．一方，Γ ′ を Γ の関数 Γ ′ (Γ ) とみなすとシンプレクティック条件を満たしながら時間発展をおΓ ′ の時間発展は以下のように計算される．ここで Γこなう分子動力学法をシンプレクティック分子動力学法という．

36

3The Molecular Simulation Society of Japan

シンプレクティック分子動力学法における時間発展

く時間発展を厳密に解くことはできないが，近似としてq と p の関数として記述される任意の物理量 A(q, p)

eDU 2 eDK ∆t eDU 2 による時間発展（あるいはより高次の鈴木・トロッター分解に基づく時間発展）は可能で∆tの時間変化は以下のように書ける．
3N (∑∂AdA=dti=1

∂qi

q̇i +

∂Aṗi
∂pi

ある．eDK ∆t および eDU ∆t はそれぞれ K および U を)

3N (∑ハミルトニアンとみなした場合の時間発展演算子であるのでこれらの時間発展演算子による時間発展はシン)

プレクティック条件 (10) 式を満たす．各部分がシンプレ∂A ∂H∂A ∂H−∂q∂p∂piii ∂qii=1[ 3N ()]∑ ∂H ∂∂H ∂=−A(q, p)∂pi ∂qi∂qi ∂pii=1=∆tによる時間発展もシンプレクティック条件を満たす．

DK , DU はその定義 (18) 式より q, p に以下のように演算する．

∂K2, DKqi = 0, DU qi = 0,∂pi
∂U2, DU
pi = 0, DK pi = 0
DU pi = −∂qi
DK qi =

1 行目から 2 行目の計算には正準方程式 (2) 式を用いた．ここで DH をDH ≡i=1

∂H ∂∂H ∂−∂pi ∂qi∂qi ∂pi)(13)

ここで DK は 2 回以上 qi に演算した場合 0 になり，DUる．これらを用いると eDK ∆t , eDU ∆t による q, p の時

dA= DH Adt(14)とあらわされる．この形式解はA(t + ∆t) = eDH ∆t A(t)(15)と書ける．eDH ∆t を時間発展演算子と呼ぶ．
ハミルトニアンが運動量のみに依存する項と座標の間変化は以下のように計算される．
()
1 2DK ∆t2eqi = 1 + DK ∆t + DK (∆t) + · · · qi2∂K∆t= qi +(21)∂pieDK ∆t pi = pi

(22)

eDU ∆t qi = qi(23)

DU ∆te

みに依存する項に分かれている場合

H(q, p) = K(p) + U (q)(16)∂Upi = pi −∆t∂qi

H=3N3N∑∑∂K ∂∂U ∂DH =−∂p∂q∂qi ∂piiii=1i=1= DK + DU(24)粒子 i の座標 r i とその運動量 pi を用いてハミルトニアンがを考える．この場合，時間発展演算子もN∑
p2i2mii=1(17)n+ 12pi= pni + F ni(18)これらの時間発展演算子 DK , DU をもちいて全体の時間発展演算子 DH を例えば∆t(19)と分解する．この分解を鈴木・トロッター分解という．
一般に分子動力学シミュレーションでは eDH ∆t∆t2に基づ(26)n+ 12pi∆tmi∆tn+ 1pn+1= pi 2 + F n+1ii2r n+1= r ni +iと定義した．

eDH ∆t = eDU 2 eDK ∆t eDU 2 + O((∆t)3 )(25)

すると r i , pi の時間発展はみに依存する項に対応する時間発展演算子をそれぞれ3N∑∂U ∂DU ≡ −∂qi ∂pii=1+ U (r)と書ける場合，(19) 式のように時間発展演算子を近似と分離できる．ここで運動量のみに依存する項と座標の∆t(20)

は 2 回以上 pi に演算した場合 0 になることに注意すと定義すると A(q, p) の時間変化は3N∑∂K ∂,DK ≡∂pi ∂qi
i=1∆t

クティック条件を満たしているので eDU 2 eDK ∆t eDU 2(12)3N (∑∆t(27)(28)のようにおこなえる．ここで上付きの n は n ステップ目での物理量をあらわし，粒子 i に働く力をFi = −∂U∂r i(29)と書いた．(26)−(28) 式の時間発展は速度ベルレ法といわれる．速度ベルレ法を変形して得られるベルレ法，

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

リープ・フロッグ法もシンプレクティック分子動力学法である．(19) 式の代わりに eDH ∆t を∆t(35) 式のように 3 つの演算子を 1 つにまとめるには
ベーカー・キャンベル・ハウスドルフの公式を 2 回繰∆teDH ∆t = eDK 2 eDU ∆t eDK 2 + O((∆t)3 )(30)と近似すると，り返して使う．
AAeC = e 2 eB e 2(39)の関係がある場合，演算子 D を

pn ∆t= r ni + imi 2n+ 1
ri 2n+ 21

pn+1= pni + F iin+ 12= ri
r n+1i(31)∆t(32)

pn+1 ∆t+ i
mi 2(33)のように時間発展できる．これは位置ベルレ法と呼ばれる [4]．速度ベルレ法と位置ベルレ法のどちらを使っても構わないが，位置ベルレ法では力を半ステップ ∆t/2ずれた時刻で計算する．このため瞬間圧力(N)N
1 ∑ p2i (t) ∑P(t) =+r i (t) · F i (t)3V
mii=1(34)i=1など座標・運動量と力の両方の関数を計算するときには同時刻のそれぞれの項（この場合は運動エネルギーとヴィリアル）をすぐには計算できない．同時刻の運動エネルギーとヴィリアルを計算する必要がある場合には二度手間だが r n+1から F n+1を計算するか，n − 12ii
ステップ目と n + 12 ステップ目のヴィリアルから内挿して近似的に n ステップ目のヴィリアルを求める．

437

シンプレクティック分子動力学法における保存量
前節では時間発展演算子 eDH ∆tAeD ≡ e 2 eBと定義する．ベーカー・キャンベル・ハウスドルフの公式を使うと

11
1
1A+B+ [A, B]+ [A, [A, B]]+ [[A, B], B]+· · ·2
44824(41)と求まる．
AeC = eD e 2(42)D=であるからもう一度ベーカー・キャンベル・ハウスド
ルフの公式を使うとC = A+B −1
1[A, [A, B]] + [[A, B], B] + · · · (43)2412となる．(35) 式と (39) 式を比較すると (43) 式より

DH̃ = DK + DU
1
(2[[DU , DK ], DK ]+24−[DU , [DU , DK ]]) (∆t)2 + · · ·1 つにまとめた演算子をDK + DU =3N∑∂H ∂i=1∆t∂qi ∂pi(45)となる．ここで 2 行目の { , } はポアソン括弧である．

DU , DK の交換子は{演算子 A, B, C の間に[DU , DK ] = −C(36)の関係が成り立つとき

11
C = A+B+ [A, B]+ {[A, [A, B]] + [[A, B], B]}+· · ·212(37)とあらわされる．ここで [A, B] は A と B の交換子で[A, B] ≡ AB − BA(38)と定義される．C は A, B およびその交換子の組み合わせであらわされる．(37) 式をベーカー・キャンベル・ハウスドルフの公式と言う．

i=1= DHキャンベル・ハウスドルフの公式を用いる．

e e =e3N∑∂H ∂(35)と定義する．DH̃ を DK , DU であらわすにはベーカー・A B∂pi ∂qi−= { , H}

∆teDH̃ ∆t ≡ eDU 2 eDK ∆t eDU 2≈ eDH ∆t(44)と求まる．また (18) 式よりを (19) 式のように近似的に分解した．分解された時間発展演算子を逆に(40),3N∑∂H ∂Hi=1∂pi ∂qi}(46)と求まる．同様に [DU , DK ] と DU , DK の交換子は
3N ∑3N ∑∂H ∂ 2 H ∂H 
[DU , [DU , DK ]] =,∂qi ∂pi ∂pj ∂qj i=1 j=1(47)3N ∑3N ∑∂H ∂ 2 H ∂H [[DU , DK ], DK ] =,
∂pi ∂qi ∂qj ∂pj i=1 j=1(48)38

The Molecular Simulation Society of Japanと計算される．(47), (48) 式を (44) 式に代入して学法を使えば精度を上げることは可能である．

• ハミルトニアンから正準方程式として導かれる運動方程式ではない場合（能勢・フーバー熱浴など）にはDH̃ = { , H}

3N3N ∑ 1∑∂H ∂ 2 H ∂H
2,+ 24∂pi ∂qi ∂qj ∂pj
i=1 j=1使えない．

3N ∑3N2∑∂H ∂ H ∂H −(∆t)2∂qi ∂pi ∂pj ∂qj i=1 j=1+ ···(49)となる．よって，(19) 式の近似に基づくこの分子動力学シミュレーションは
3N 3N1  ∑ ∑ ∂H ∂ 2 H ∂H
H̃ ≡ H +224∂pi ∂qi ∂qj ∂pj

仮想時間 t′ における能勢の熱浴 [5, 6] の運動方程式（本連載第 1 回の (9)−(11) 式 [9]）は能勢のハミルトニ
アン HN から導かれる正準方程式である．そのためハプレクティック分子動力学法を導くことは可能である．
しかし，実時間 t での時間発展のほうが便利なことが3N ∑3N∑∂H ∂ 2 H ∂H −(∆t)2∂q∂p∂p∂qiijji=1 j=1(50)をハミルトニアンとみなしたときの時間発展を厳密に解く方法であるといえる．H̃ を影のハミルトニアンといい，これは厳密に保存する．このため長時間シミュレーションをおこなっても誤差が蓄積しにくい．一方，
シンプレクティックではない方法を用いた場合，H̃ のような保存量は存在しない．このため長時間シミュレーションをおこなうとハミルトニアンがだんだん異なる多い．実時間 t における能勢の熱浴のシンプレクティッ
ク形式が能勢・ポアンカレ熱浴である [2, 3]．
能勢・ポアンカレ熱浴のハミルトニアンは{N∑ p′2iHNP = s+ U (r)22m
isi=1}
P2+ s + gkB T0 log s − H02Q

= s {HN (r, p′ , s, Ps ) − H0 }(52)で与えられる．ここで H0 は HN の初期値である．Psは能勢の熱浴の自由度 s に対する正準共役な運動量である．Q は熱浴の質量に相当し s の時間発展を制御す値になってしまう．
また (19) 式や (30) 式では時間発展演算子を時間反転に対して対称に分割している．このためこのシンプレクテック分子動力学法は時間反転可逆な時間発展アルる．kB はボルツマン定数，T0 は設定温度である．g は系の自由度をあらわすパラメーターで N 個の原子からなる系では g = 3N とする．このハミルトニアンから運動方程式は以下のように導かれる．

ゴリズムにもなっている．すなわち

p′i
mi s
ṗ′ i = sF iPsṡ = sQN∑ p′2iṖs =− gkB T0ms2ii=1ṙ i =exp(−DH̃ ∆t) exp(DH̃ ∆t)= exp(DH̃ ∆t) exp(−DH̃ ∆t)= 15 温度一定のシンプレクティック分子動力学法：能勢・ポアンカレ熱浴

ミルトニアン HN を分割して仮想時間 t′ におけるシンi=1 j=1+···• アルゴリズムの開発が必ずしも簡単ではない．

(51)が成り立っている．
シンプレクティック分子動力学法の特徴をまとめると以下のようになる．
長所• 予測子・修正子法に比べると計算精度があまり良くない．通常用いられる方法ではハミルトニアンの誤差（グローバルエラー）は (∆t)2 に比例する．ただし，労力はかかるがより高次のシンプレクティック分子動力(54)
(55)
(56)

HN は保存されるので (56) 式で• 簡単で安定性が良い．
• ハミルトニアンの保存が良い．
• 時間反転に対し対称である．
短所

(53)

HN − H0 = 0(57)という関係を用いた．
能勢・ポアンカレ熱浴における時間発展は以下のように導かれる．ハミルトニアン HNP を以下のように 3つに分割する．

HNP = HNP1 + HNP2 + HNP3(58)分子シミュレーション研究会会誌 “アンサンブル” Vol. 11, No. 1, January 2009 (通巻 45 号)(N)∑ p′2i
HNP1 = s+ gkB T0 log s − H0 (59)
2mi s2i=1
HNP2 = sU (r)P2
HNP3 = s s2Q(74) 式では
1= 1 − x + x2 − x3 + · · ·1+x(60)(61)

各ハミルトニアンよる時間発展は{N ()∑ ∂HNP1 ∂∂HNP1 ∂·DHNP1 r i =−·∂p′i
∂r i
∂r i
∂p′ii=1}∂HNP1 ∂
∂HNP1 ∂+−ri∂Ps ∂s∂s ∂Psp′i(62)=mi sN∑p′2i
DHNP1 Ps =− gkB T0 log s + H0 − gkB T0
2mi s2i=139(76)という関係式を用いた．各演算子による時間発展がわかったので (68) 式の分割の順番で (69)−(74) 式の関係を用いることにより能勢・ポアンカレ熱浴の時間発展をおこなうことができる．ここで左向きの矢印 ← は計算プログラムでの代入を意味する．
∆teDHNP3 2 の演算：()2Ps ∆ts ← s 1+2Q 2/()Ps ∆tPs ← Ps
1+2Q 2(77)
(78)

∆t

(63)
DHNP2 p′i = sF i

(64)

DHNP2 Ps = −U (r)Ps
DHNP3 s = sQ
P2
DHNP3 Ps = − s2Q

(65)

∆t2
∆t
Ps ← Ps − U (r)2p′i ← p′i + sF i(66)

(79)
(80)

eDHNP1 ∆t の演算：
(67)
ri ← ri +と計算される．時間発展演算子を例えば∆t∆teDHNP ∆t = eDHNP3 2 eDHNP2 2 eDHNP1 ∆t× eeDHNP2 2 の演算：

DHNP2 ∆t
DHNP3 ∆t22ePs ← Ps +3+ O((∆t) ) (68)p′i∆tmi s( N∑ p′2i

2mi s2i=1(81)− gkB T0 log s)+H0 − gkB T0 ∆tのように分割し (62)−(67) 式の関係を用いると各変数(82)の時間変化は以下のように求まる [3]．

p′er i = r i + i ∆t
(69)
mi s( N∑ p′2ieDHNP1 ∆t Ps = Ps +− gkB T0 log s22misi=1)
DHNP1 ∆t+H0 − gkB T0 ∆teDHNP2 ∆t p′i = p′i + sF i ∆t(70)(71)

Ps = Ps − U (r)∆t(72)2PPseDHNP3 ∆t s = s + s ∆t + s s2 (∆t)2Q
4Q)2(Ps∆t= s 1+(73)2Q
P2P3eDHNP3 ∆t Ps = Ps − s ∆t + s2 (∆t)2 + · · ·2Q
4Q
/()
Ps
∆t= Ps1+(74)
2Qe

DHNP2 ∆t(73) 式では(DHNP3sPs22Q∆teDHNP2 2 の演算：∆t2∆tPs ← Ps − U (r)2p′i ← p′i + sF i

(83)
(84)

∆teDHNP3 2 の演算：()2Ps ∆ts ← s 1+2Q 2/()
Ps ∆t
Ps ← Ps
1+2Q 2(85)(86)

6 能勢・ポアンカレ熱浴と能勢・フーバー熱浴の比較：
水分子への応用
前回の解説記事で能勢・フーバー熱浴 [5-7] による時間発展 [8] を述べた（本連載第 1 回の (49)−(54) 式 [9]）．
そこで用いた変数 r i , pi , ζ は正準変数ではない．また)

その時間発展はシンプレクティック条件 (10) 式を満た=0(75)さない．一方，能勢・ポアンカレ熱浴で用いる変数 r i ,40

The Molecular Simulation Society of Japan
0.15

7 まとめ

δH(t) (kcal/mol)(a) Nose-Poincare thermostat

今回はシンプレクティック分子動力学法について解0.10

説した．シンプレクティック分子動力学法では影のハミルトニアンという厳密に保存する量が存在する．そ0.05のため安定にシミュレーションをおこなえる．実時間における能勢の熱浴のシンプレクティック形式である

0.00

能勢・ポアンカレ熱浴についても説明した．能勢・フーバー熱浴を用いた場合，ハミルトニアンの値は時間と-0.05ともに増大してしまう．これに対し，能勢・ポアンカ

δH(t) (kcal/mol)(b) Nose-Hoover thermostat

レ熱浴を用いた場合，ハミルトニアンの値は良く保存0.10する．そのため能勢・フーバー熱浴を用いるよりも能勢・ポアンカレ熱浴を用いる方が時間ステップ幅を長0.05く設定できる．能勢・ポアンカレ熱浴に限らずシンプ
レクティック分子動力学法は今後ますます使われてい0.00-0.050.0くであろうと考えている．

0.5

1.01.5

参考文献[1] H. Yoshida, Phys. Lett. A, 150, 262 (1990).

t (ns)

図 1: ハミルトニアンの初期値からのずれ δH(t) の時間変化．(a) 能勢・ポアンカレ熱浴，(b) 能勢・フーバー
熱浴．

p′i , s, Ps は正準変数であり，その時間発展 (77)−(86) 式はシンプレクティック条件を満たしている．この節ではこれら 2 つの熱浴を用いた場合にそれぞれハミルトニアンがどの程度保存されるかを示す [10]．
周期境界条件を課した立方体セルの中に剛体水分子を 80 個用意し，温度を 300 K，密度を 0.997 g/cm3 に設定した．静電相互作用の計算にはエワルド法を用い[2] S. D. Bond, B. J. Leimkuhler, and B. B. Laird, J. Comput. Phys., 151, 114 (1999).
[3] S. Nosé, J. Phys. Soc. Jpn., 70, 75 (2001).
[4] M. Tuckerman, B. J. Berne, and G. J. Martyna, J.
Chem. Phys., 97, 1990 (1992).
[5] S. Nosé, Mol. Phys., 52, 255 (1984).
[6] S. Nosé, J. Chem. Phys., 81, 511 (1984).
[7] W. G. Hoover, Phys. Rev. A, 31, 1695 (1985).
[8] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M.
L. Klein, Mol. Phys., 87, 1117 (1996).

ミュレーションをおこなった．剛体分子の時間発展に[9] 奥村久士，分子シミュレーション研究会会誌 “アンサンブル”，10, No. 4, 29 (2008).

はシンプレクテック解法を用いた [11]．

[10] H. Okumura, S. G. Itoh, and Y. Okamoto, J. Chem.

た．時間ステップ幅を ∆t = 4 fs として t = 1.5 ns シ

能勢の熱浴のハミルトニアンの初期値からのずれ

δH(t) = HN (t) − H0 .

(87)の時間変化を図 1 に示す．能勢・ポアンカレ熱浴を用

Phys., 126, 084103 (2007).
[11] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, D. Newns, and G. J. Martyna, J. Chem. Phys.,116, 8649 (2002).

いた場合では影のハミルトニアンが存在するので，時奥村久士 (博士 (理学))：〔経歴〕

間ステップ幅を ∆t = 4 fs にとってもハミルトニアン
著者紹介

し能勢・フーバー熱浴ではそのような保存量が存在し2002 年慶應義塾大学理工学研究科博士課程修了，東京大学日本ないので時間がたつにつれてハミルトニアンはその初学術振興会特別研究員 (PD)，分期値からずれてしまう．そのため能勢・フーバー熱浴

子科学研究所助手，名古屋大学を使う場合には時間ステップ幅を ∆t = 4 fs には取れ

COE 特任講師を経て現在ラトガース大学研究助教授．
〔専門〕
統計力学，分子シミュレーショがその初期値から大きくずれることはない．それに対ず，∆t = 0.5 − 2 fs 程度に設定することが多い．

ン．
〔趣味〕海外旅行，落語鑑賞．