26

The Molecular Simulation Society of Japan

分子動力学シミュレーションにおける温度・圧力制御
第 5 回：パリネロ・ラーマンの方法，圧力一定のガウス束縛法，
圧力一定のベレンゼンの方法
自然科学研究機構
分子科学研究所

計算科学研究センター

理論・計算分子科学研究領域

奥村久士

1

はじめに

hokumura@ims.ac.jpz

前回は圧力を制御する手法としてアンダーセンの方y法を解説した．この方法は立方体のシミュレーション
セルを等方的に膨張・縮小させることにより圧力を制御する．今回は圧力を制御するもう 1 つの拡張系の方法であるパリネロ・ラーマンの方法を解説する．このc手法はシミュレーションセルの変形を許す手法である．
拡張系の方法以外にも圧力を制御する方法はある．
本連載第 3 回において温度を制御するガウス束縛法とbベレンゼンの方法について解説した．今回は圧力を制御するガウス束縛法とベレンゼンの方法についても解ax説する．

2

パリネロ・ラーマンの方法
パリネロ・ラーマンの方法 [1] では図 1 に示すように平行六面体セルを用いる．この平行六面体セルの各辺を a, b, c とする．各辺のベクトルを横に並べた行列を
ax bx cx()L = a b c =  ay by cy (1)azbzczションセル.と書ける．ここで任意の行列 A と B について

AT B T = (BA)(6)

スケールした座標 r̃ と行列 L を用いてラグランジ

r i = Lr̃ i(2)アン

L=とスケールし，速度をṙ i = Lr̃˙ i(3)cTa·c

アンという．ここで第 3 項の Tr (· · ·) は行列の対角和である．この項は( T ) 1
1
M Tr L̇ L̇ = M (a2 + b2 + c2 )22b·c c·cN( T )11∑Tmi r̃˙ i Gr̃˙ i −U (r̃, L)+ M Tr L̇ L̇ −P0 V
2 i=12(7)を導入する．これをパリネロ・ラーマンのラグランジとスケールする．ここで L を用いて計量テンソル
 T a() T T
G=L L =  b  a b ca·a a·b
=  b·a b·b
c·a c·b

Tの関係が成り立つことを用いた．と書く．L を用いて座標を

図 1: パリネロ・ラーマンの方法におけるシミュレー(8)と書き換えることができ，シミュレーションセルの運(4)を導入する．計量テンソル G を用いると速度の 2 乗は()TTT
˙Lr̃˙ i = r̃˙ i LT Lr̃˙ i = r̃˙ i Gr̃˙ i (5)ṙ 2i = ṙ Ti ṙ i = Lr̃ i

動エネルギーを意味する．このラグランジアンから r̃ iと L に関する運動方程式を導出する．
ラグランジュの運動方程式はd ∂L∂L=˙dt ∂ r̃ i∂r̃ i(9)分子シミュレーション研究会会誌 “アンサンブル” Vol. 11, No. 4, October 2009 (通巻 48 号)

d ∂L∂L=dt ∂ L̇
∂L(10)と書ける．ここで粒子 i に働く力をFi = −とあらわされる．以下，これらの式の各項を計算する．
まず G は対称行列であることに注意する．

( )T)T= LT L = G
GT = LT L = LT LT(11)

この点に注意すると，ラグランジアン (7) 式の第 1 項にあらわれる粒子の速度の 2 乗は()Tr̃˙ i Gr̃˙ i =r̃˙ i1 r̃˙ i2 r̃˙ i3
G11 G12 G13r̃˙ i1×  G12 G22 G23   r̃˙ i2 G23G33∂U∂r i(21)と書いた．

(G1327

次にシミュレーションセルの運動方程式を導出する．

(10) 式左辺は以下のように計算できる．
( T )d ∂L
1 d ∂= MTr L̇ L̇
dt ∂ L̇
2 dt ∂ L̇
1 d ∂ ( 2= Mȧ + ȧ2y + ȧ2z + ḃ2x + ḃ2y + ḃ2z
2 dt ∂ L̇ x)+ċ2x + ċ2y + ċ2z
= Mr̃˙ i3222= G11 r̃˙ i1 + G22 r̃˙ i2 + G33 r̃˙ i3+ 2G12 r̃˙ i1 r̃˙ i2 + 2G13 r̃˙ i1 r̃˙ i3 + 2G23 r̃˙ i2 r̃˙ i3(12)dL̇ = M L̈dt(22)

これに対し (10) 式右辺は{ N}
∂L
∂
1∑T
˙
˙
=mi r̃ i Gr̃ i − U (r̃, L) − P0 V
∂L
∂L 2 i=1(23)と書ける．以下，(23) 式右辺の各項を計算する．(23)と書ける．

()∂ (˙T ˙ )r̃ i Gr̃ i = 2 G11 r̃˙ i1 + G12 r̃˙ i2 + G13 r̃˙ i3
˙∂ r̃i1()= 2 Gr̃˙ i(13)1Tであるから r̃˙ i Gr̃˙ i の r̃˙ 微分は∂ (˙T ˙ )r̃ i Gr̃ i = 2Gr̃˙ i∂ r̃˙ i(14)となる．ここで L = (a b c) の逆行列を用いて行列 σと計算できる．このことから (9) 式左辺はd (を)d ∂L=mi Gr̃˙ i = mi Ġr̃˙ i + mi Gr̃¨ idt ∂ r̃˙ idt(15)と求まる．さらに (9) 式右辺は∂L
∂U (r̃, L)∂r T ∂U∂U=−=− i= −LT∂r̃ i∂r̃ i∂r̃ i ∂r i
∂r i(16)と計算できる．(9), (13), (16) 式より r̃ i の運動方程式は

1∂Ur̃¨ i = − G−1 LT− G−1 Ġr̃˙ imi∂r i(17)と書ける．ここで( )−1G−1 = L−1 LT

式右辺第 1 項は (1), (4) 式より∂ (˙T ˙ )∂ (222r̃ i Gr̃ i =a·a r̃˙ i1 + b·b r̃˙ i2 + c·c r̃˙ i3∂L∂L+ 2a·b r̃˙ i1·r̃˙ i2 + 2a·c r̃˙ i1·r̃˙ i3)+ 2b·c r̃˙ i2·r̃˙ i3() T(24)= 2 Lr̃˙ i r̃˙ iL−1 =(b × c)T

1 1 T (c × a)T  ≡ σVV(a × b)Tと定義する．すなわち(σ = b×cc×aa×b(25))(26)

である．σ を用いると (24) 式は( N)N()()T1∑
1 ∑∂T˙˙mi r̃ i Gr̃ i =mi Lr̃˙ i Lr̃˙ i σ∂L 2 i=1V i=1(27)と書ける．

(18)

次にポテンシャルエネルギー U (r̃, L) の L 微分を考える．ここではテンソル計算がやや複雑なので r̃ i の αであるから

G−1 LT = L−1(19)なので粒子の運動方程式は

1 −1r̃¨i =L F i − G−1 Ġr̃˙ i
mi成分を r̃iα と書き L の (α, β) 成分を Lαβ と書くことにする．α, β, · · · は 1–3 の値をとる．この記号法を用いると座標は(20)

riα =3∑β=1Lαβ r̃iβ(28)28The Molecular Simulation Society of Japanと書ける．ここで 1 は 3 × 3 の単位行列である．(20),と書ける．さらに U (r̃, L) の L 微分は(36) 式を数値的に時間発展させることによりシミュレーションセルの変形を許す圧力一定の分子動力学シミュN ∑3∑∂U (r̃, L)∂riγ ∂U (r)
−= −
∂Lαβ
∂L
αβ ∂riγ
i=1 γ=1=3N ∑∑レーションをおこなうことができる．また能勢の方法とパリネロ・ラーマンの方法を組み合わせることによδαγ r̃iβ Fiγ

りシミュレーションセルの変形を許す温度・圧力一定のi=1 γ=1=N∑

分子動力学シミュレーションをおこなうこともできる．

Fiα r̃iβ(29)3 温度・圧力一定のガウス束縛法

i=1−1となる．Lを用いると

r̃iβ =

本連載第 3 回 [2] で温度を制御するガウス束縛法を3∑( −1 )Lrβγ iγ(30)γ=1−∂U (r̃, L)=∂Lαβi=1 γ=1御することができる [3, 4]．ここではアンダーセンの方法と同じように立方体シミュレーションセルを等方的と書けるのでこれを (29) 式に代入すると
N ∑3∑

解説した．ガウス束縛法で温度だけでなく，圧力も制に膨張・圧縮させる．ミクロカノニカル分子動力学シ()
Fiα L−1 βγ riγ

ミュレーションの運動方程式pi
miṗi = F iṙ i =N ∑3∑

1Fiα riγ σγβ=V i=1 γ=1(31)(37)(38)に η, ζ を未定定数として以下のようにつけ加える．

すなわち

∂U
1−=∂LV(N∑)
F i rTipi+ ηr imiṗi = F i − ζpi − ηpiṙ i =σ(32)i=1となる．ここで括弧 (· · ·) 内はテンソル積であることに注意する．

(39)
(40)

瞬間温度，瞬間圧力が一定になるように η, ζ を決める．
ここで η の意味を考えておく．アンダーセンの方法

次に体積 V の L 微分を考える．体積 V は1と同様に L = V 3 を用いて座標 r i ，運動量 pi をV = a · (b × c)

1= ax by cz − ax bz cy + ay bz cx− ay bx cz + az bx cy − az by cx

r i = V 3 r̃ i(41)− 31(42)pi = Vp̃i(33)とスケールする．(41) 式の時間微分をとるとであるのでその L 微分は
by cz − bz cy az cy − ay cz ay bz − az by∂V=  bz cx − bx cz ax cz − az cx az bx − ax bz ∂Lbx cy − by cx ay cx − ax cy ax by − ay bx()=(34)
b×c c×a a×b =σ1ṙ i = V 3 r̃˙ i +V̇pL̇
ri = i + ri
3VmiL

(43)

ここでアンダーセンの方法において
1pi = mi ṙ i = mi V 3 r̃˙ i(44)と書けることを用いた．(39) 式と (43) 式を比較すると，と書ける．ゆえに∂(−P0 V ) = −P0 σ∂Lη は V の運動を決めていることがわかる．
(35)
V̇ = 3ηV(45)となる．

(10), (22), (27), (32), (35) 式よりシミュレーションセ
ルの運動方程式は[ {N}
N()()T ∑1 1 ∑L̈ =mi Lr̃˙ i Lr̃˙ i +F i rTi
M V i=1i=1{z}|]−P0 1 σ

圧力テンソル(36)

またL̇ri(46)
Lであり，ηr i は系の膨張，収縮にともなう流れにのって
ηr i =粒子が移動することを示している．

ζ を決めるには瞬間温度
T (t) =N1 ∑ p2i
gkB i=1 mi(47)分子シミュレーション研究会会誌 “アンサンブル” Vol. 11, No. 4, October 2009 (通巻 48 号)

が一定値 T0 になるように束縛をかける：

これらをまとめるとガウス束縛法による温度・圧力一定の分子動力学シミュレーションの運動方程式はN∑

p2imii=1

gkB T0 =(48)

ここで g は系の自由度であり，N 個の粒子からなる系· ṗi = 0i

mii=1(49)(40) 式を代入して
N∑pV̇ = 3ηVN
∑p(57)
(58)
(59)

iN
∑

· Fimiζ = i=1N
∑ p2−η =−pi ∂U·mi ∂r ii=1
gkB T0i· (F i − ζpi − ηpi ) = 0i

mii=1pi+ ηr i
miṗi = F i − ζpi − ηpiṙ i =の場合には g = 3N −1 にする．この両辺を微分すると，N∑pi=1(50)

N∑pN
∑pi

η = − N∑· Fimiζ = i=1N∑ p2−η(51)−ηimi·N ∑N∑pj∂U∂2U+·ri ·∂r i i=1 j=1∂r i ∂r j mjN ∑N∑∂U∂2U+ri ·· r j + 9P0 V∂r i i=1 j=1∂r i ∂r j

ri ·i=1(61)ii=1miとなる．

を得る．

この運動方程式により系を時間発展させた場合の分η を決めるには瞬間圧力(N)N∑
1 ∑ p2i∂UP(t) =−ri ·3V i=1 mi i=1∂r i

布関数 f (r, p, V ) は(52)

3P0 V =N∑

p2ii=1mi−N∑ri ·i=1∂U∂r iii=1imi)gi− kB T02m2ii=1{z}運動量に束縛(53)両辺の微分をとるとN∑
2p · ṗf (r, p, V ) = δ(N∑ p2|が一定値 P0 になるように束縛をかける．

3P0 V̇ =(60)mii=1

よって，29( ())
N∑∂U
1
gkB T0 −ri ·− P0 V× δ3∂r ii=1|{z}座標に束縛NN N∑∂U ∑ ∑∂2U−ṙ i ·−ri ··ṙ j∂r i i=1 j=1 ∂r i ∂r ji=1{}
U (r) + P0 V
× exp −
kB T0(62)(54)(49) 式より右辺第 1 項は 0 である．さらに (39), (45) 式であることが知られている [5]．定温定圧アンサンブルを使って変形すると，分布関数が得られるが，これは定温定圧アンサンブル

9ηP0 V = −N∑pi=1imi·ではない．なぜならば運動量への束縛だけでなく座標N∑∂U∂U−ηri ·∂r i∂r ii=1にも束縛がかかっているからである（(62) 式右辺第 2項）．ガウス束縛法で温度制御をおこなうと，運動量N ∑N∑pj∂2U·−ri ·∂r∂rmijji=1 j=1−ηN ∑N∑
i=1 j=1ri ·空間ではカノニカルアンサンブルを得られなくても座標空間ではカノニカルアンサンブルを得られた．しかし，ガウス束縛法で圧力制御までおこなうと座標空間

2∂ U· rj∂r i ∂r jと同じように exp [{−U (r) + P0 V } /kB T0 ] に比例する(55)ですら定温定圧アンサンブルを得ることはできない．

4 圧力一定のベレンゼンの方法

よって，N∑
pi=1

η=− N∑

imi·N ∑N
∑pj
∂U
∂2U+ri ··∂r i i=1 j=1∂r i ∂r j mj

圧力を制御するベレンゼンの方法もある [6]．この方法では瞬間圧力N ∑N∑∂2U∂U+ri ·· r j + 9P0 Vri ·∂r i i=1 j=1∂r i ∂r ji=1(56)1P(t) =3V(N∑ p2N∑

∂U−ri ·mi i=1∂r ii=1i)
(63)

30

The Molecular Simulation Society of Japan

5 まとめが設定圧力 P0 にdP
P(t) − P0=−(64)dtτと近づくように座標と体積をスケールする．ここで τは P が P0 に近づく時間スケールを決定する緩和時間である．アンダーセンの方法や圧力を制御するガウス束縛法と同じように η を用いて系の膨張，収縮速度をあらわし，今回は圧力を制御する手法としてパリネロ・ラーマンの方法を解説した．この方法はアンダーセンの方法に比べると複雑ではあるが，結晶の相転移や構造変化など結晶格子の形状が変化する場合によく用いられている．
さらに温度と圧力を制御するガウス束縛法とベレンゼンの方法も説明した．通常はこれら 2 つの方法より

pṙ i = i + ηr imiV̇ = 3ηV(65)

も能勢・アンダーセンの方法を使うのが一般的である．
ガウス束縛法やベレンゼンの方法で温度と圧力を制御(66)しなければいけないということはまずないが，古いプと時間変化させる．η を決定するために圧力の時間変ログラムの中にはこれらの手法を使っているものもあ化を等温圧縮率る．その場合は正しい定温定圧アンサンブルを得るこ

β=−1V(∂V∂P)とができない点に注意して使わなければならない．

(67)T

参考文献

を用いてあらわすと，瞬間圧力の時間変化はdP
1 dV
3η
=−=−dtβV dtβ(68)と書ける．これを (64) 式と比較するとβ(69)(P − P0 )3τとなる．これを (65) 式に代入すると圧力を制御するベ
η=レンゼンの方法の運動方程式pβṙ i = i +(P − P0 )r imi3τ

[1] M. Parrinello and A. Rahman, Phys. Rev. Lett. 45(1980) 1196.
[2] 奥村久士，分子シミュレーション研究会会誌 “アンサンブル”，11, No. 2, 43 (2009).
[3] D. J. Evans and G. P. Morriss, Chem. Phys. 77 (1983)
63.
[4] D. J. Evans and G. P. Morriss, Phys. Lett. 98A (1983)

(70)

433.
[5] S. Nosé, J. Chem. Phys., 81, 511 (1984).

が得られる．ここで等温圧縮率 β は正確にわかっている必要はない．なぜならば β の不確定性は緩和時間 τに組み込めるからである．
温度を制御するベレンゼンの方法においてシミュレーションの途中で運動量をスケールしたように（本連載[6] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak, J. Chem. Phys., 81,3684 (1984).
[7] T. Morishita, J. Chem. Phys., 113, 2976 (2000).

第 3 回 (3) 式 [2]），この手法においても以下のように座標と体積をスケーリングすることにより (70) 式を近奥村久士 (博士 (理学))：〔経著者紹介

似的に実現できる：

r i → µr i

(71)

V → µ3 V(72)歴〕2002 年慶應義塾大学理工学研究科博士課程修了，東京大学日本学術振興会特別研究員 (PD)，分子科学研究所助手，名古屋大学 COE 特任講

ここで µ はµ=1+または{µ=β∆t(P − P0 )3τ} 13β∆t1+(P − P0 )τ

師，ラトガース大学研究助教(73)

所准教授〔専門〕統計力学，(74)

であたえられる．
ベレンゼンの方法で温度と圧力を制御しても定温定圧アンサンブルを得ることはできない．詳しくは森下の論文を参照してほしい [7]．

授を経て現在，分子科学研究
分子シミュレーション．
〔趣味〕
海外旅行，落語鑑賞．