22

The Molecular Simulation Society of Japan

分子動力学シミュレーションにおける温度・圧力制御
第 4 回：アンダーセンの方法と能勢・アンダーセンの方法
自然科学研究機構
分子科学研究所

計算科学研究センター

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

奥村久士

1

ࡒࠢࡠࠞࡁ࠾ࠞ࡞(N, V, E)

はじめにhokumura@ims.ac.jp࿶ജ৻ቯ(N, P, H)

これまで分子動力学シミュレーションで温度を制御する方法を解説した．今回は圧力を制御するアンダーᢿᾲო

センの方法を解説する．また温度と圧力を両方制御する方法である能勢・アンダーセンの方法も解説する．

2アンダーセンの方法ቅ┙♽ࡇࠬ࠻ࡦߦ⚿߮ઃߌࠄࠇߚ♽物理系のハミルトニアンは

H0 =N
∑
p2
i

2mi
i=1+ U (r)(1)とあらわされる．ここで i 番目の粒子の座標，運動量，質量をそれぞれ r i , pi , mi とあらわし，N 個の粒子の座標の集合を r ≡ {r 1 , r 2 , · · · , r N } と書いた．
体積を大きくすれば圧力は下がり，体積を小さくすれば圧力は上がる．よって圧力を制御するには図 1 のようにピストンをつけて体積 V を変化させればよい．
ただし図 1 では横方向にのみ伸縮させているが実際の

図 1: 圧力一定の方法の模式図.の関係がある．
なぜ (3) 式で第 2 項を無視して良いのだろうか？理由の 1 つは運動エネルギーのうち温度に寄与するのは乱雑な熱運動のみだからである．互いに近くにいる粒子の速度の平均で決まる局所的な流れの寄与は運動エネル
1

ギーから取り除かなくてはいけない．(3) 式を L = V 3を使って書きかえるとṙ i = Lr̃˙ i + L̇r̃ i|{z} |{z}分子動力学シミュレーションでは 3 次元的に等方的に体積を増減させる．このように体積を変化させて圧力を制御する分子動力学法は H. C. Andersen によって提案された [1]．アンダーセンの方法では立方体セルの一1辺の長さ L = V 3 を用いて座標を
1r i = V 3 r̃ i(2)1V̇ṙ i = V r̃˙ i +V 3 r̃ i3Vとなる．この第 2 項を無視し
1ṙ i = V 3 r̃˙ i(3)PV 2+U(Vr̃)++ P0 V
HA =22M
2mi V 3i=1
1
3(4)(5)はスケールされた運動量で pi と1への寄与が異なってしまう．このようなセル中の位置(3) 式の第 2 項を無視する理由である．ただし，一度ṙ i = Lr̃˙ i としてハミルトニアン HA を導入した後，正準方程式を導く際には ṙ i = Lr̃˙ i + L̇r̃ i として計算しなければならない．
ハミルトニアン HA より運動方程式を導くと，を導入する．M はピストンの仮想的な質量である．p̃ipi = p̃i /V 32 項は一様な体積の増減にともなう流れをあらわしている．このため，第 2 項は無視すべきなのである．また，シミュレーションセル中の粒子はどこにあってもに依存する項が存在することは望ましくない．これもミルトニアンp̃2iとなる．ここで第 1 項は乱雑な熱運動をあらわし，第を取り入れてしまうと，セル中の位置 r̃ i によって速度ピストンを以下のように付けることにより拡張系のハN∑流れ同等に扱われなければいけない．もし (3) 式で第 2 項とスケールする．この時間微分をとると
13熱運動

(7)

(6)

∂HAdr̃ ip̃i
==2dt
∂ p̃i
mi V 3dp̃i∂HA∂U=−=−dt∂r̃ i∂r̃ idV∂HAPV==dt∂PVM(8)(9)(10)分子シミュレーション研究会会誌 “アンサンブル” Vol. 11, No. 3, July 2009 (通巻 47 号)

dPV∂HA=−dt∂V(N)N∑1 ∑ p̃2i=+F i · r i − P0
3V i=1 mi V 23i=1

3 能勢・アンダーセンの方法
温度と圧力を制御して定温定圧アンサンブルを得るには能勢の方法 [2, 3] とアンダーセンの方法を組み合わせれば良い．物理系のハミルトニアン (1) 式に熱浴(11)となる．ここで F i は i 番目の粒子に働く力でFi ≡ −∂U∂r iとピストンの自由度をつけ加えて拡張系のハミルトニアン(12)HNA =と定義される．これを r̃ i , V に関する方程式に書き直N∑p̃′2i122mi s2 V 3i=1|{z+ U (V 3 r̃)}

物理系すと，

2V̇ ˙r̃ i(13)1 Fi −
3V
mi V 3{(N)}N∑
2 2
11 ∑˙3V̈ =mi V r̃ i +F i · r i −P0
M 3V i=1i=1{z}|r̃¨ i =231

瞬間圧力 P(t)(14)+を得る．ここでは時間と空間のスケールを両方おこなっている．時間のスケールについては ′ をつけ，空間のスケールについては˜をつけてあらわしている．これらの変数の関係は
1r i = V 3 r̃ i

圧力 P(t) をあらわしている．すなわち，

pi = p′i /s = p̃i /V(15)∆t = ∆t /s

圧力 P0 より低ければ，体積の加速度 V̈ が負になりである．

(V̈ < 0)，系の収縮速度を大きくさせて瞬間圧力 P(t)を上げる．逆に P(t) が P0 より高ければ，体積の加速方程式を求めるとて P(t) を下げる．このフィードバック機構により瞬間圧力 P(t) の値を設定圧力 P0 の周りに揺らがせる．
ピストンの質量 M の値を大きくすると V はゆっくり時間変化し，逆に M の値を小さくすると V は速く時間変化するようになる．瞬間圧力 P の時間変化と同じくらいのスケールで V の値も時間変化するように Mの値を調節すると効率よく圧力を制御できる．
運動方程式 (8)−(11) 式あるいは (13), (14) による時間発展においてハミルトニアン HA は保存量である．

HA = 一定(16)平衡状態では体積は平均値の周りでゆっくり揺らぐだ
けで，体積の運動エネルギー PV 2 /2M はハミルトニア(20)13′と書きかえられる．これは系の瞬間圧力 P(t) が設定度 V̈ が正になり (V̈ > 0)，系の膨張速度を大きくさせ

(21)
(22)

∂HNAp̃′i(23)=2′∂ p̃i
mi s2 V 3dp̃′i∂HNA∂U=−=−(24)dt′∂r̃ i∂r̃ ids∂HNAPs==(25)′dt∂PsQ(N)dPs∂HNA
1 ∑ p̃′2i=−=− gkB T0 (26)dt′∂ss i=1 mi s2 V 23
dr̃ i=dt′dV∂HNA
PV==(27)dt′
∂PVM
∂HNA
dPV
=−′dt
∂V(N)
N
∑
1 ∑ p̃′2i=+F i · r i − P0
3V i=1 mi s2 V 32i=1(28)PV 2(17)≪ HA2M
このためエンタルピー H は近似的に一定値とみなせる．
(18)

よってアンダーセンの方法では近似的に N P H 一定のアンサンブル（等エンタルピー等圧アンサンブル）が得られる．

= p̃′i /sV13

正準方程式から r̃ i , p̃′i , s, Ps , V , PV についての運動

ン HA に比べて小さな値をとることが多い．

H ≡ H0 + P0 V ≈ 一定

(19)

ピストン

熱浴となる．(14) 式 { } 内の第 1，2 項の和は系内部の瞬間

M V̈ = P(t) − P0 .

Ps2
PV 2+ gkB T0 log s ++ P0 V
2Q
2M {z}|{z} |となる．さらにこれを実時間 t で時間発展するように書きかえる．座標，運動量をそれぞれ，r̃ i から r i へ

p̃′i から pi へ変数変換してdr ipV̇= i +ridtmi3V()dpiṡ
V̇= Fi −+pi
dts 3V(29)
(30)

24

The Molecular Simulation Society of JapandsPs= sdtQN∑dPsp2i=− gkB T0dt
mii=1(31)4 定温定圧アンサンブルが実現することの証明
ハミルトニアン HNA は一定値 E のままであり，拡(32)dV
PV= s(33)
dtM)}
{(NN∑
dPV
1 ∑ p2i= s+
F i · r i − P0dt
3V i=1 mi i=1(34)を 得 る ．こ れ ら の 運 動 方 程 式 (23)−(28) 式 ま た は(29)−(34) 式にしたがって時間発展させることにより定温定圧アンサンブルでの分子動力学シミュレーショ

張系全体ではミクロカノニカルアンサンブルになる．
そのため拡張系全体の分配関数は∫ ∞ ∫ ∞∫ ∞∫ ∞∫∫Y =dsdPsdVdPV dr̃ dp̃′0dr̃dp̃′ = s3N drdp

ここでも粒子の運動量と体積にフィードバックがかる．物理系の瞬間温度−∞0と書ける．ここで仮変数から実変数への変数変換 r i =11V 3 r̃ i , pi = p̃′i /V 3 s をおこなう．r̃ i , p̃′i の微分はンを行うことができる.
かってそれぞれ温度と圧力を制御していることがわか

−∞

( ()
1p̃′P23× δ H0 V r̃,+ s + gkB T0 log s12QsV 3)
2PV+ P0 V − E+(39)2M(40)と書ける．ここで x, y, z の各成分について

drix = V 3 dr̃ix , dpix = dp̃′ix /sV 3(41)driy = V dr̃iy , dpiy = dp̃′iy /sV(42)

1N
∑p2i

1T (t) =gkB i=1 miを使って (31)，(32) 式を書きかえると
( )
d ṡ
gkB
=
(T (t) − T0 )
dt s
Q(35)

(36)

が設定温度 T0 より高ければ ṡ/s は大きくなって pi を小さくし，温度を下げようとする．逆に，瞬間温度 T (t)が設定温度 T0 より低ければ ṡ/s は小さくなって pi を大きくし，温度を上げようとする．ここで抵抗係数とは言っても ṡ/s は T (t) に依存して正にも負にもなることに注意する．

(37)

i=1を使って (34) 式を書きかえると(38)となる．すなわち瞬間圧力 P(t) が設定圧力 P0 より高ければ体積の運動量 PV は大きくなって体積 V を大きくし，圧力を下げようとする．逆に，瞬間圧力 P(t) が設定圧力 P0 より低ければ体積の運動量 PV は小さくなって体積 V を小さくし，圧力を上げようとする．このように粒子の運動量と体積にフィードバックをかけることにより温度と圧力を制御している．

1
3

1dpiz = dp̃′iz /sV 3(43)る．この変数変換をおこなうと∫ ∞ ∫ ∞∫ ∞∫ ∞∫∫Y =dsdPsdVdPV dr dp0−∞−∞0(P23N× s δ H0 (r, p) + s + gkB T0 log s2Q)
2
PV++ P0 V − E2M(44)となる．ここで δ 関数の恒等式

δ(f (s)) =dPV
= s {P(t) − P0 }dt

13なので (40) 式右辺の係数が s3N であることに注意すとなる．ここで ṡ/s は (30) 式右辺第 2 項，運動量に対i=113driz = V dr̃iz ,する抵抗係数に対応している．物理系の瞬間温度 T (t)また物理系の瞬間圧力(N)N∑
1 ∑ p2iP(t) =+F i · ri
3Vmi1

δ(s − s0 )|f ′ (s)|ただし

f (s0 ) = 0(45)を用いる．今の場合

Ps2PV2

−− P0 V  E − H0 (r, p) −
2Q 2Ms0 = exp
gkB T0(46)であるから，δ 関数の積分を実行すると∫ ∞ ∫ ∞∫ ∞∫ ∞∫∫Y =dsdPsdVdPV dr dp0−∞0−∞/()
gkB T0× s3Ns

PV2Ps2
−− P0 V E − H0 (r, p) −

2Q 2M

× δ s−exp

gkB T0

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

∫ ∞

∫ ∞
∫ ∞
∫
∫
1dPsdVdPV dr dp
gkT0
B−∞0−∞{(23N + 1P× exp
E − H0 (r, p) − s
gkB T02Q)}
PV 2×−− P0 V
2M=

HNPA1 = s×−∞

∫ ∞dPs(47)(× exp

定温定圧アンサンブル

dPV1gkB T0−∞)
E − Ps2 /2Q − PV 2 /2M

kB T0(N∑p̃′2ii=12mi s2 V 32

HNPA2 = s

HNPA4 = s(54)Ps22Q(55)

各ハミルトニアンによる時間発展は{ N ()∑ ∂HNPA1 ∂∂HNPA1 ∂·DHNPA1 r̃ i =−·∂r̃ i∂r̃ i∂ p̃′i∂ p̃′ii=1+(48)=には本連載第 1 回 [4] と同様に，g = 3N とすると定温DHNPA1 Ps =定圧アンサンブルが得られる．
能勢・アンダーセンの方法の時間発展(29)−(34) 式の運動方程式を差分化し実際に分子動力学シミュレーションをおこなう手法として，ここでは能勢・ポアンカレ熱浴 [5-7] を用いた方法を紹介する．
能勢・ポアンカレ熱浴は実時間における能勢の熱浴のシンプレクティック解法を与える．能勢・ポアンカレ
熱浴は本連載第 2 回に詳しく解説しているので，もし興味があればそちらも参照してほしい [8]．能勢・ポア
ンカレ熱浴とアンダーセンの方法を組み合わせたハミ

1+ U (V 3 r̃)}P2
P2+ s + gkB T0 log s + V + P0 V − H02Q
2M{}
= s HNA (r̃, p̃′ , s, Ps , V, PV ) − H0(49)で与えられる．ここで H0 は HNPA の初期値である．こ∂HNPA1 ∂∂HNPA1 ∂−∂Ps ∂s∂s ∂Ps∂HNPA1 ∂∂HNPA1 ∂+−∂PV ∂V∂V∂PV得られることがわかる．実時間 t でサンプルする場合

p̃′2i2 22msi V 3i=1(53)2M{}
1HNPA3 = s U (V 3 r̃) + P0 V

合には g = 3N + 1 とすれば定温定圧アンサンブルがルトニアンは{N∑
HNPA = s+ gkB T0 log s − H0PV2と計算される．このことから仮想時間 t′ で平均する場5)(52)となる．ここで，g = 3N + 1 とすると{}∫∫∫ ∞H0 (r, p) + P0 V
dV dr dp exp −Y =kB T0|0{z}∫ ∞25p̃′i}r̃ i(56)

2mi sV 3N∑p̃′2i

22mi s2 V 3i=1− gkB T0 log s+ H0 − gkB T0N∑p̃′2iDHNPA1 PV =2 53msi V 3i=1(57)
(58)

PV22M
PV
DHNPA2 V = sM
1DHNPA3 p̃′i = sV 3 F i

DHNPA2 Ps = −(59)
(60)
(61)
13

DHNPA3 Ps = −{U (V r̃) + P0 V }()N1 ∑F i · r i − P0
DHNPA3 PV = s
3V i=1(62)
(63)

PsQ
P2
DHNPA4 Ps = − s2Q
DHNPA4 s = s(64)(65)と計算される．時間発展演算子を例えば∆t∆teD∆t = eDHNPA4 2 eDHNPA3 2のハミルトニアン HNPA から正準方程式を求めると実× eDHNPA2 2 eDHNPA1 ∆t eDHNPA2 2

時間における運動方程式 (29)−(34) 式が導かれる．た× eDHNPA3 2 eDHNPA4 2 + O((∆t)3 ) (66)∆t∆t

だし，その際に∆t∆tのように分割し，この順番で (56)−(65) 式の関係を用

HNPA − H0 = 0(50)の関係を用いる．
ハミルトニアン HNPA による時間発展を導くためにHNPA を以下のように分割する．
HNPA = HNPA1 + HNPA2 + HNPA3 + HNPA4 (51)いると定温定圧アンサンブルにおける時間発展を以下のように得ることができる．ここで左向きの矢印 ← は計算プログラムでの代入を意味する．
∆teDHNPA4 2 の演算：()2Ps ∆ts ← s 1+2Q 2(67)26The Molecular Simulation Society of Japan/(Ps ← Ps1+Ps ∆t2Q 2)(68)

保つよりも圧力を一定に保っておこなわれることが多1∆tp̃′i ← p̃′i + sV 3 F i(69)2][1∆tPs ← Ps − U (V 3 r̃) + P0 V(70)2()
N1 ∑∆tPV ← PV + s(71)F i · r i − P0
3V i=12eの演算：

PV2∆t
2M 2PV ∆tV ← V +sM 2e(72)Ps ← Ps +い実験も多い．実験条件とシミュレーション条件を直接比較して議論するにはアンダーセンの方法を用いて圧力を制御するのが便利である．

参考文献

p̃′imi sV( N∑23∆t(74)p̃′2i[2] S. Nosé, Mol. Phys., 52, 255 (1984).
[3] S. Nosé, J. Chem. Phys., 81, 511 (1984).

(73)の演算：r̃ i ← r̃ i +い．特に近年は実験手法も進歩し圧力を加えた興味深[1] H. C. Andersen, J. Chem. Phys. 72 (1980) 2384.

Ps ← Ps −DHNPA1 ∆t

さらにこれらの分子動力学シミュレーションをおこなうための具体的な手順も述べた．実験は体積を一定に∆teDHNPA3 2 の演算：

DHNPA2 ∆t2する方法として能勢・アンダーセンの方法も解説した．

[4] 奥村久士，分子シミュレーション研究会会誌 “アンサンブル”，10, No. 4, 29 (2008).
[5] S. D. Bond, B. J. Leimkuhler, and B. B. Laird, J. Comput. Phys., 151, 114 (1999).

− gkB T0 log s (75)[6] S. Nosé, J. Phys. Soc. Jpn., 70, 75 (2001).

)[7] H. Okumura, S. G. Itoh, and Y. Okamoto, J. Chem.
Phys., 126, 084103 (2007).

22mi s2 V 3i=1+H0 − gkB T0 ∆t(76)[8] 奥村久士，分子シミュレーション研究会会誌 “アンPV ← PV +N∑
i=1p̃′2i
3mi sV

53∆t(77)∆teDHNPA2 2 の演算：
P 2 ∆tPs ← Ps − V2M 2PV ∆tV ← V +sM 2

著者紹介(78)(79)

歴〕2002 年慶應義塾大学理工学研究科博士課程修了，東究員 (PD)，分子科学研究所助手，名古屋大学 COE 特任講∆t師，ラトガース大学研究助教

1∆tp̃i ← p̃i + sV 3 F i(80)2{} ∆t1Ps ← Ps − U (V 3 r̃) + P0 V(81)2()
N1 ∑∆tPV ← PV + s(82)F i · r i − P0
3V i=12∆t

eDHNPA4 2 の演算：

6奥村久士 (博士 (理学))：〔経京大学日本学術振興会特別研eDHNPA3 2 の演算：()2Ps ∆ts ← s 1+2Q 2)/(Ps ∆tPs ← Ps
1+2Q 2サンブル”，11, No. 1, 35 (2009).

(83)(84)

まとめ
今回は圧力を制御する分子動力学法としてアンダーセンの方法を解説した．また温度と圧力の両方を制御

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