統計力学（小正準集合の方法）
能勢修一
慶應義塾大学理工学部・大学院理工学研究科

1
分子動力学（MD) 法は相互作用をしている多数の粒子の集団の示す性質を調べる統計力学的シミュレーション法の一つであり、粒子の運動を古典力学の運動方程式を数値的に解くことにより求めている。この点では、万有引力の下で運動する太陽系の惑星の動きを調べることや打ち上げられたロケットの運動を調べることと同じことである。しかし、その背景となる考え方はかなり異なっている。惑星やロケットの運動では、個々の質点の動きを長期間にわたって正確に求めることが必要であり、初期条件によりどのように質点の運動が変わるかが重要なテーマとなる。
これに対し、統計力学的なシミュレーションでは、粒子の集団が初期状態でどのようであったかを忘れ、初期状態に依存しない平衡状態に達することがまず必要となる。その上で、多くの粒子配置、粒子の運動の統計的な平均として粒子の集団の示す巨視的な性質を求めていく。個々の粒子の運動を非常に正確に求めることよりも統計平均して得られた値の信頼度を高くすることの方が重要となる。
シミュレーションである物理量（例えば圧力）を知るためにはどのような量の平均を求めればよいのかは、統計力学が教えてくれる。統計力学は熱力学で取り扱うような巨視的な体系の性質を、原子分子の間に働く相互作用および力学法則という微視的な情報より求める手段を与えてくれるものである。温度が一定の状態で考えるか、全エネルギー一定で考えるかなど、調べようとする体系の置かれている条件により、具体的に統計力学を適用するときの手順が異なっていくる。どのような物理量の組（一種類の粒子のみからなる系
では通常三つでよい）を一定に保って実験またはシミュレーションを行うかによって用いるべき統計力学的集合（アンサンブル）が指定される。古典力学に基づく MD シミュレーションでは（保存系であれば）体系の全エネルギーが保存される。これに対応する統計力学の集合はミクロカノニカル（小正準）集合である。この集合では全エネルギー E 、体積V 、粒子数 N が一定の条件に対応するカノニカル（正準）集合と、温度 T 、体積 V 、化学ポテンシャル µ が一定の条件に対応するグランドカノニカル（大正準）集合がよく用いられる。
多くの場合、小正準集合に比べ、正準集合の方が取り扱いが容易なため、多くの統計力学の教科書では小正準集合については基本原理を述べただけで、具体的な物理量の表現は正準集合の場合のみ記している。MD シミュレーションを行う場合には、これでは不便なので小正準集合について少し詳しく説明する。
物理量を統計力学に従い求める時注意すべきことは、求めるべき量によりすべての統計力学的集合で同じ表式となる場合と、集合ごとにその表式が異なる場合があることである。体系が熱平衡状態にあるかどうかを判定する基準となる最も基本的な熱力学量を熱力学ポテンシャルと呼ぶ。結論を先に述べておくと、この熱力学ポテンシャルの変数および2一次微分量についてはすべての統計力学的集合で同じ表式となるが、２次以上の微分量に関係した量ではその表式は集合により異なっている。
小正準集合で考えている体系は他の系と相互作用をしない孤立系であり、体系の全エネルギーは一定に保たれている。孤立系ではエントロピーは増大するのみというよく知られた熱力学第２法則の表現（エントロピー増大則）より熱平衡ではエントロピーが最大となる。従って、全エネルギー E 、体積 V が一定の条件（小正準集合に対応）の熱力学ポテン
シャルはエントロピー S である。温度 T 、体積 V が一定の場合（正準集合に対応）の熱力学ポテンシャルはヘルムホルツの自由エネルギー F = E − T S であり、非平衡状態ではF は減少していく。また、温度 T 、圧力 P 一定の場合の熱力学ポテンシャルはギブスの自由エネルギー（この量を単に熱力学ポテンシャルと呼ぶ時もある）G = E − T S + P Vである。
統計力学の基本原理は、微視的な状態についての和の形で表される分解関数を巨視的な熱力学量の代表である熱力学ポテンシャルと関係づけることである。有名なボルツマンの関係式 S = k ln W （小正準集合）やその変形である F = −kT ln Z （正準集合）がこれである。ここで、W (E, V, N ) は全エネルギーが E となる微視的状態の数であり、Z(T, V, N )は正準集合の分配関数 Z = i exp[−Ei /kT ] である。このように微視的な量と巨視的な量の対応関係が確立すると、その他の巨視的な熱力学量に対する微視的な統計力学的表式は熱力学ポテンシャルの全微分を基に熱力学関係式を用いることにより導くことができる。
熱力学ポテンシャルの変数となる量（状態変数）および一次微分量については、すべての集合で統計力学的な表現が同じになる。熱力学ポテンシャルの全微分の式に現れる量がこれにあたる。正準集合では dF = −SdT − P dV + µdN より、状態変数が T 、V 、N 、一次微分量が S 、P 、µ である。全エネルギー E も E = F + T S の関係より一次微分量までで表すことができる。
しかし、２次微分量（これは物理量の揺らぎに関係している）およびさらに高次の微分量の表式は集合により異なってくる。注意すべき点は密度を一定のまま粒子数を無限大とする熱力学極限においては、これらの熱力学量の統計力学的表式は集合により異なっているがその値は同じになることである。代表的な例は、定積熱容量 CV = (∂E/∂T )V である。正準集合では、エネルギーの揺らぎと結び付いたよく知られた表式CV =(δE)2 N kT 2(正準集合)となる。この表式は正準集合でのみ用いることのできるものである。もし、この表式をそのまま小正準集合に適用するとエネルギー保存よりエネルギーの揺らぎはないため定積熱容量は 0 であるという結果になる。しかし、これは正しくない。小正準集合では次の式

3を用いるべきである（ただし、全エネルギーが運動エネルギー K =2i pi /2mi の部分と位置エネルギー Φ の部分に分離できる。また体系が単原子系であるとしている）。

CV =k1 − (1 − 2/3N )K1/K

(小正準集合)粒子数 N が大きい熱力学極限では

CV =k2/3N − (δK)2 /K2と運動エネルギー K （または δK = −δΦ より位置エネルギー Φ についても同じ結果を与える）の揺らぎと関係づけられる。
古典統計力学での小正準集合

N 個の粒子からなる古典的な体系を考える。この時、統計力学で考える一つの微視的状態はすべての粒子の運動量（まとめて p = (p1 , p2 , · · · , pN ) と表す）、と座標（同様にq = (q1 , q2 , · · · , qN ) と表す）により定義される位相空間 Γ = (p, q) の一点で表される。時間とともに各粒子がそれぞれ動いて粒子配置が変っていくことは、Γ の空間での点の移動として表される。この点が動いてできる曲線を位相空間での軌道を呼ぶ。エネルギー保存則を満たさねばならないので軌道は Γ 空間のどの点でも通ることができるわけではない。
許される軌道は Γ 空間の中の H(p, q) = E を満たす超平面の上に制限されている。古典統計力学ではある物理量 A の平衡状態での値は位相空間で定義されるある関数 B(p, q) の軌
道に沿っての平均（これを時間平均と呼ぶ）

1 t0 +T
A = B̄ =B(t)dt
T t0により与えられるとする。通常は、T → ∞ の極限をとり始めの時刻 t0 に依存しない場合を考える。
もし、Γ 空間内の軌道が H(p, q) = E を満たす超平面のほとんどすべての点を通るようになると上の時間平均は H(p, q) = E を満たすすべての状態にわたっての平均B = B(p, q)ρ(p, q)dpdq
B(p, q)δ(H((p, q) − E)dpdq=δ(H((p, q) − E)dpdqに等しくなる（これを位相平均と呼ぶ）。上の式に現れる ρ(p, q) は位相空間 Γ 上での確率分布を表す分布関数である。小正準集合で許されるエネルギーの幅を無限小とした極限の4形が δ 関数を用いた表現

ρ(p, q) = Cδ(H(p, q) = E)である。
小正準集合の方法はこの位相平均に基づいて定式化されている。本当に時間平均 B̄ と位相平均 B が等しくなるかどうかは難しい問題である。この等式が成り立つとするのがエルゴードの仮説である。またこの等式が成り立つとき体系はエルゴード的であるという。自由粒子や調和振動子の集まりのように運動方程式が解析的に積分可能な場合にはエルゴード性は成り立たない（自由度の数と同じ数だけの独立な保存量が存在するとき運動方程式は可積分であるという）。多粒子系で相互作用がある程度複雑な場合にはエルゴード的になっているものと仮定して我々はシミュレーションを行うが、場合によっては注意が必要である（低温の結晶系ではソリトン的な動きが観測されることがある）。
位相空間での確率分布の形を知るには古典力学のリュービルの定理を利用する。
位相空間で点 Γ の近傍 dΓ（qi が qi と qi + dqi 、pi が pi と pi + dpi の間の値を持つ領域）に運動の状態が存在する確率が f (Γ)dΓ となるように分配関数 f (Γ) を定義する。
微小体積 dΓ については確率の保存則より∂f+ div (Γ̇f (Γ)) = 0∂t∂fdΓdt は、その時間の間に微小体∂t
積の表面を通して変化する量 −div (Γ̇f (Γ))dΓdt に等しいことをいっている。ここでが成り立つ。これは、微小体積内での dt の間の変化

div (Γ̇f (Γ)) = ∇ · (Γ̇f (Γ)) = (∇ · Γ̇)f (Γ) + Γ̇ · ∇f (Γ)
を用いると∂f+ Γ̇ · ∇f (Γ) = −(∇ · Γ̇)f (Γ)∂tと変形できる。正準方程式で表されるような通常の力学では ∇ · Γ̇ = 0 となる（保存系の時）。

  ∂ ṗi∂ q̇i∇ · Γ̇ =+∂pi ∂qii  ∂  ∂H ∂H∂=−+=0∂pi∂qi∂qi ∂pii

これより

∂fdf=+ Γ̇ · ∇f (Γ) = 0dt∂t5となる。

df∂fはストークス微分とも全微分とも呼ばれる量である。 dΓdt∂tが位相空間の一点（Γ）での時間変化を表すのに対して、全微分は位相空間の中の軌道にここで定義した時間微分沿って動いた微小体積の間での変化を見ている。
df= 0 は古典力学に従う運動では、軌道に沿っての状態分布の形はリュービルの定理
dt変化しないことを示している。t = t0 に位相空間の点 Γ0 に f (Γ0 )dΓ の密度で状態があると、t = t1 には Γ1 に状態は動いており、そこでの分布密度 f (Γ1 )dΓ は t = t0 の場合と同じである。

f (Γ0 ) = f (Γ1 )(または f (Γ(t0 )) = f (Γ(t1 )) )ここで、もしエルゴード性が満たされており、位相空間の中の許される状態のほとんどすべてを通るとすると、位相空間の中の許される状態では、おなじ体積には等しい確率密度で状態が分布していることになる。

• δ 関数の分布になることを書く。

小正準集合でまず計算すべき量は全エネルギーが E と E + ∆E の間の微視的状態の数である。これは、∆ → 0 の極限では



ω(E) =δ(H(Γ) − E)dΓである。これを熱力学のエントロピー S とS(E, V, N ) = k ln ω(E)と結び付けるのがボルツマンの関係式である。ここでは、よく用いられる W (E) の代わりにエントロピーの２つの定義を区別して表すために ω(E) または Ω(E) の記号を用いる。
ボルツマンの関係式は、微視的に求められる量 ω(E) と対象の巨視的性質を表す熱力学量（今の場合はエントロピー S ）との対応を与える統計力学の最も基本的な関係である。ただし、状態の数については、エネルギー E 以下のすべての微視的状態の数 E
Ω(E) =ω(E )dE = θ(E − H(Γ))dΓを用いたほうがいろいろの表式が簡単でわかりやすくなる。θ(x) は Heavyside の階段関数

θ(x) =

1x>00x<06でこの関数の微分がデルタ関数となる (θ (x) = δ(x))。これからは、主にこちらの表式を用い、エントロピーはS  = k ln Ω(E)と表されると考える。エネルギーの増加とともに微視的状態の数が非常に急激に増大する
通常の場合には、ω(E) を Ω(E) に置き換えたために生ずる違いは体系に含まれる粒子の数を N とすると、1/N 程度であり、粒子数 N → ∞ となる極限ではその差はほとんど影響しない。ただし、小正準集合でのいろいろの量の統計平均を計算するには、エネルギーを E の値に限った分布関数 δ(H(Γ) − E) を用いる必要がある。
A = A(Γ)δ(H(Γ) − E) dΓ
これからは、古典粒子系で、全エネルギーがH(p, q) =

 p2ii2mi+ Φ(q)と書ける場合に限ることにする。この時、小正準集合の分配関数および、統計平均量を比較的簡単に導くことができる。

P = (p1 , p2 , · · · , pN ) という運動量の空間を考え、X =数変換する。この時の体積要素は dX = √ EΩ(E) =dE
 EdE=
 E=

dE 

p1p2pN√ , √ ,··· , √m1m2mNdPである。
Πi mi

δ(H(Γ) − E  ) dΓdQ

dP δiidQ

 p2dX2mi

Πi mi δ(+ Φ(q) − E 

Xi2 + Φ(q) − E  )i

ここで x = |X| を考えると
 EΩ(E) =dE
dQ dxx3N −1 A3N δ(x2 + Φ(q) − E  ) E=dE
 E

dQ



dxx3N −1 A3N

δ(x −2E  − Φ(q))
E  − Φ(q)11dQA3N (E  − Φ(q)) 2 (3N −2)23N
1= dQA3N(E − Φ(q)) 23N=dE 

と変7と位置座標について運動エネルギー K(q) = E − Φ(q) の 3N/2 乗の積分という形で表すこ√3N −1とができた。この式に現れる定数 A3N は Πi mi と 3N 次元の球の表面積をと3N x
C∞
3N/22π表すときの数因子 C3N =の積である（Γ(x) はガンマ関数 Γ(x) =y x−1 e−y dyΓ(3N/2)0である）。
 1Ω(E) = A3N(E − Φ(q))3N/2 dQ
3N

ω(E) では、この表式が1dΩ(E)
ω(E) == A3NdE2 (E − Φ(q))3N/2−1 dQとなる。
種々の熱力学量についての小正準集合での統計力学的な表式は、この分配関数 Ω(E) を微分することにより、導くことができる。小正準集合では、エネルギー E 、体積 V 、および粒子数 N が状態の様子を指定する変数となる。エントロピー SA3N(E − Φ(q))3N/2 dq
S(E, V, N ) = k ln Ω(E) = k ln
3N Vを熱力学の全微分式

dS =1PµdE + dV − dNTTTに従って微分していくことにより、種々の熱力学量の微視的な表現が求められる。例えば、温度の逆数 1/T はエントロピー S をエネルギー E で微分することにより求められる。
3N
3N/2−1(E − Φ(q))
dqk(E − Φ(q))3N/2−1 dq
3N k1
∂S2

===T
∂E2
3N/2dq(E − Φ(q))(E − Φ(q))3N/2 dq=

13N k2 E − Φ(q)mc

ここで、E − Φ(q)mc は小正準集合での運動エネルギーの平均値であるが座標 q の関数として表されるので K(q)mc のように表すこともある。ここで得られた結果は、運動エネ
ルギー K(q) = E − Φ(q) の平均が 3N/2kT に等しい

K(q)mc = E − Φ(q)mc =という等分配則を示している。

3NkT2

8
圧力はエントロピーを体積で微分した

∂S
P==kT
∂V


∂
(E − Φ(q))3N/2 dq∂V (E − Φ(q))3N/2より導かれる。
この表式のままでは積分範囲に体積の値が含まれており、微分を行うのが困難となるので、次のような手法を用いる。体積の変化は全系を一様に拡大または縮小することにより表されるとする。ここでは、考えている体系の形は立方体としておく。粒子の座標を立方体の辺の長さでスケーリングして表すと

q = V 1/3 sとなる。体積要素は dq = V N ds と変換され、積分範囲から体積の大きさが消える。

3N/2(E − Φ(q))
dq = (E − Φ(V 1/3 s))3N/2 V N dsV1

この表式で体積による微分を実行すると


∂
∂
3N/2N
1/33N/2(E − Φ(q))
dq =(E − Φ(V s))
dsV∂V V∂V1
N −1= NV(E − Φ(V 1/3 s))3N/2 ds1 ∂Φ qi3N(E − Φ(V 1/3 s))3N/2−1 −ds+VN∂qi 3V1 2i
ここでもう一度変数 q に戻すと

N=V
V(E − Φ(q))3N/2 dq

1 3N−3V 2 ∂ΦV∂qii· qi (E − Φ(q))3N/2−1 dqとなる。これより ∂Φ ∂Φ
· qi  i ∂qkNk N  i ∂qi · qi kNPi=−=1−TV2 V E − Φ(q)
V2 · 3NkT2となる。従って、圧力 P の表式は

1P =V



1N kT −3 ∂Φi∂qi

· qi9または

 p21P =3Vimii−

 ∂Φi∂qi

· qiという、ビリアル定理による表式となる。熱力学ポテンシャル（この場合は S ）の一階微分である温度および圧力の統計力学的な表式は集合に依存しないので正準集合の場合と同じ形である。
ゆらぎと関係した２階微分量になると集合による違いが出てくる。
度 E で微分すると

∂2S
∂E 2

11 ∂T=− 2=− 2T ∂ETV∂E∂T∂S
1= をもう一∂ETVとなる。定積熱容量 CV はエネルギーの温度微分で定義される。これをづけると

CV =

∂E∂T

1=− 2TV



∂2S
∂E 2
V

 2 ∂P∂ Sは ∂E∂Vと関係づけられる。
∂T V V 2 ∂ S∂P1 ∂PP ∂T
∂
=− 2=
∂E∂V
∂E ∂T VT ∂E VT∂E V∂PP
1−=∂T V CV T
CV T 2∂P
1=−PT
CV T 2∂T Vとなる。同様にして圧力係数

グリュナイゼン係数 γγ=Vを用いると上の関係式はとも表される。 また、圧縮率∂P∂E∂2S∂E∂VV=CVV

=∂P∂TVPγ−TVCV T 2 2 ∂P∂ Sはと関係づけられる。
∂V V∂V 2 E 2  ∂ SP∂
1 ∂PP ∂T==−∂V 2 E ∂V T E T ∂V E T 2 ∂V E

∂2S
∂E 2に関係V10
ただし、エネルギー一定の条件での圧縮率は、実験値との比較が難しいので、これを断熱圧縮率で表すと
∂P∂V


=E∂P∂V+SPT CV∂P∂T
=V∂P∂V+SPγTVとなる。
定積熱容量の場合について具体的に計算を示すと次のようになる。


3N k
3N/2−1 2 dq (E − Φ(q))
∂ S∂
2=
∂E 2 V∂E (E − Φ(q))3N/2 dq (E − Φ(q))3N/2−2 dq
3N k 3N=−1 
22(E − Φ(q))3N/2 dq2
3N/2−1
(E − Φ(q))
dq 
3N k 3N
−2 2 (E − Φ(q))3N/2 dq  
213N113N k 3N−1=−k22K mc Kmc2
K2mc
これより、定積熱容量の小正準集合での表式は  −13N kT 12CV = k 1 − 1 −3N2K mc  −112=k 1− 1−
K
3NK mcとなる。通常用いる表式は、粒子数 N が大きく、運動エネルギー K の揺らぎ δK が Kmcに比べて小さい場合のものである。運動エネルギー K の逆数の平均は
 
11=K mc
Kmc + δK mc


11=
δK
Kmc 1 + K
mcmc


2δK1
δK=+− ···
1−Kmc
Kmc
Kmc

1δK 2 mc+ ···=1+Kmc
K2mc11と表されるので、これを CV の表式に代入すると
−12δK 2 mc
CV = k 1 − 1 −1+3N
K2mc−12δK 2 mc−=k3N
K2mcと運動エネルギー K の揺らぎと関係づけた表式となる。これは、最初 Lebowitz 等によっ



∂P
1
グリュナイゼンの係数 γ および断熱圧縮率の逆数 BS == −Vについての
κS∂V S
小正準集合の表式は   ∂Φ∂ΦNk
3N
1−1 V−γ=−VV2∂V∂V mc K mc
 2 N kT∂ Φ +VBS = Nk∂V 2 mc
V 1 + 2γ −CV  2 3N∂Φ∂Φ
1−−1 V−2∂V∂V mc Kて求められたよく知られた表式である。

mcとなる。そして、N → ∞ の時の表式はそれぞれ、 NkV∂Φ
γ=+ 3δδKCV∂VN kT 2mc2  2 2 N kTV∂ Φ∂Φ +VBS = −
δ2Nk
∂V mc kT∂Vmc
V 1 + 2γ −CVとなる。以上の式で1  ∂Φ∂Φ  ∂Φ ∂qi==·· qi∂V∂q∂V3V∂qiiiiの略記号を用いた。同様に∂2Φ∂=2∂V∂V∂Φ∂V1  ∂Φ· qi
3V i ∂qi 2 1∂ 2Φ
1  ∂Φ
· qi +qj ·· qi=− 2
3V i ∂qi3V∂qj ∂qii,j=である。

∂∂V