はじめに
(古典)分子動力学シミュレーションでは原子(イオン，分子)間に働く力をモデル化し多粒子系のニュートン方程式を解き，粒子群の時間発展を記述します。力の情報が分かると系の圧力をビリアル定理から直接計算することができます。特に，力が二体力ポテンシャル(二つの粒子の位置のみで決まる保存力を与える粒子缶相互作用エネルギー)で与えられる場合には系の圧力を容易に計算することができ，
それによって系の状態方程式を決定することもできます。
本稿では圧力の計算方法について説明し，それに基づいてプロトタイプコードをPythonによって作成します。
用いる2体ポテンシャルはよく知られたレナード-ジョーンズポテンシャルです。
これらは分子動力学シミュレーションやモンテカルロ法のコードの一部に組み込むことができます。
初歩的な統計力学の知識のみ仮定しました。
当該分野をはじめて学ぶ読者の方の参考になれば幸いです。
＃数式・内容に間違いがあればご指摘いただければありがたいです。。
どんな人に役立ちそう？
● 分子動力学シミュレーションやモンテカルロ法における圧力の計算法の初歩について知りたい方
● 状態方程式の計算法について理論的な興味がある人
状態方程式
相互作用ポテンシャルが一定値の場合，原子間に力が働かず，系は理想気体の状態方程式
$$p =\frac{Nk_BT}{V} {\tag 1}$$
に従います。
原子間相互作用がある場合には下記のファンデルワールスの状態方程式
$$(p+\frac{N^2 a}{V^2})(V-Nb) =Nk_BT \
\iff p =\frac{Nk_BT}{V-b} -\frac{N^2 a}{V^2} {\tag 2}$$
が有名です。$a$および$b$がゼロのときに理想気体の状態方程式に一致します。
$b$は原子(分子)の大きさを考慮し，体系の体積を補正する補正するためのパラメータです。
一方，$a$は原子間に相互作用による圧力の実効的な増加補正のためのパラメータです。
本稿では原子間ポテンシャルを用いて第一原理的に$a$の表式を導くことも行います。
圧力の導出
本稿のメインです。
ここでは統計力学により，粒子数($N$)・温度($T$)・体積($V$)一定下における圧力$P(N,V,T)$(状態方程式)の微視的な表式を導きます。
系のヘルムホルツの自由エネルギーを$F(N,V,T)$とし，分配関数を$Z(N,V,T)$とするとすると，熱平衡状態における圧力$P$は，
$$P= -(\frac{\partial F(N,V,T)}{\partial V})_{N,T}$$
$$= k_B T (\frac{\partial ln Z(N,V,T)}{\partial V})_{N,T}{\tag 3}$$
とかけます。
ここで，$F$と$Z$の間に成り立つ関係式
$$F = -k_B T ln \ Z(N,V,T) {\tag 4}$$
を用いました。
$Z(N,V,T)$は，量子統計力学では
$Z(N,V,T) =\sum_{n} \exp({ -E_n/(k_B T)}){\tag 5}$
と表されます。ここで$E_n$は量子状態$n$における量子力学的エネルギー期待値です。
古典力学では，半量子論的取り扱い(古典統計力学近似)により
$${{ Z(N,V,T) = \frac{1}{N!(2 \pi \hbar)^{3N}} \int\int.... \int e^{-H(p,q)/k_BT} }} {\mathrm d^3} p_1{\mathrm d^3} p_2...{\mathrm d^3} p_N \, {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N $$
$$ = \Lambda \int e^{-U(q)/k_BT} {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N {\tag 6}$$
となります。ここで，分母の$N!$は同種粒子が区別できない性質(不可弁別性)を考慮するための補正項です。$H$は系のハミルトニアンであり，(位相)空間における粒子$i$の(正準)運動量を$\bf{p_i}$，(正準)座標を$\bf{q_i}$として粒子間相互作用ポテンシャルエネルギーを$U(V,\bf{q_1,q_2,..., q_N})$と表すと，$H= \sum_{i}\frac{|p_i|^2}{2m_i}+U(\bf{q_1,q_2,..., q_N})$となります。
$\Lambda$は$V$に依存しない量です。
式(6)を式(3)に代入し，最終的に$P$の統計力学的な表式が得られます。
$$P=\frac{Nk_BT}{V}-\frac{\int e^{-U/k_BT}(\frac{\partial U}{\partial V})\ {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N}{\int e^{-U/k_BT}\ {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N} {\tag 7}$$
右辺第二項は，原子間ポテンシャル$U$の体積$V$による変化率の熱平均値に比例する物理量となっていることが分かります
絶対零度のときの圧力が$P(T=0K)=\frac{\partial U}{\partial V}$であることを思い出すと，この式が原子間相互作用による圧力の発生を示していることがイメージしやすいと思います。
ここまでは，同種粒子系を仮定したことを除き，古典統計力学的に厳密な話です。
吟味
もし$U$が一定ならどうなるでしょうか。その場合，粒子間に相互作用がないため，理想気体の状態方程式(1)に帰着するべきです。実際にそうなるのか，(7)式から確かめてみましょう。このとき$\frac{\partial U}{\partial V}=0$となるので，式(7)の右辺第二項の分子はゼロになり，結局，$P=\frac{Nk_BT}{V}$となって期待したとおり理想気体の状態方程式(1)となることがわかります。
また，式(7)とファンデルワールスの状態方程式(2)とを比べることで，$b<<1$のときのファンデルワールス状態方程式中のパラメータ$a$の統計力学的表現として
$$a = \frac{V^2}{N^2}\frac{\int e^{-U/k_BT}(\frac{\partial U}{\partial V})\ {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N}{\int e^{-U/k_BT}\ {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N}{\tag 8}$$
を得ます。これにより，$a$が確かに原子間ポテンシャルの情報を含んでいることが分かります。
これ以上の解析は$U$の具体的な関数形を知らないと分かりませんが，$U$が二体ポテンシャルの場合は表式が著しく簡単化され，数値計算を行うことも難しくありません。
2体ポテンシャルに対する圧力の表式
原子$i$と原子$j$の原子間相互作用が両原子の原子間距離$r_{ij}$のみに依存する場合は，系の全相互作用ポテンシャルエネルギーは，原子のペアが感じる相互作用ポテンシャル(ペアポテンシャル)$U_{ij}=U_{pair}(r_{ij})$の和になります。
$$U=\sum_{i(\ne j)}U(r_{ij}){\tag 9}$$
このとき，原子$j$が原子$i$に及ぼす力$f_{ij}$は，
$$f_{ij} = -\frac{\partial U_{pair}(r_{ij})}{\partial r_{ij}}\frac{\bf r_{ij}}{r_{ij}} {\tag
{10}}$$
となるので，式(7)の右辺第二項の分子の$\frac{\partial U}{\partial V}$は，
$$\frac{\partial U}{\partial V}=\sum_{i(\ne j)}\frac{\partial U_{pair}(r_{ij})}{\partial V} =\sum_{i(\ne j)}\frac{\partial U_{pair}(r_{ij})}{\partial r_{ij}} \frac{r_{ij}}{3V}\
=\frac{1}{3V} \sum_{i(\ne j)}\frac{\partial U_{pair}(r_{ij})}{\partial r_{ij}} r_{ij}\
= -\frac{1}{3V} \sum_{i(\ne j)} r_{ij} \cdot f_{ij}{\tag {11}}$$
となるので，2体ポテンシャルによる圧力の表式は，熱平均値を表す記号$\langle \cdot\cdot\cdot \rangle$を用いて
$$P=\frac{Nk_BT}{V}-\frac{\int e^{-U/k_BT}(-\frac{1}{3V} \sum_{i(\ne j)} r_{ij} \cdot f_{ij})\ {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N}{\int e^{-U/k_BT}\ {\mathrm d^3} q_1 {\mathrm d^3} q_2 ...{\mathrm d^3} q_N} \
= \frac{Nk_BT}{V}+\frac{1}{3V}\langle \ \sum_{i(\ne j)} \bf r_{ij} \cdot f_{ij} \ \rangle {\tag {12}}$$
となることが分かります。
この$\sum_{i(\ne j)} \bf r_{ij} \cdot f_{ij} \ \rangle $がビリアルとよばれる量です。
この式の第二項とファンデルワールスの状態方程式(式2)を比較すれば，$a$の意味が相互作用によるものであることがはっきりと分かります。
Pythonによる数値計算
ここでは2体ポテンシャルとして中性原子系を良く記述することがしられている"レナード・ジョーンズポテンシャル"をとり，実際にビリアルを計算するプロトタイプコードを作成します。
原子間距離を$r_{ij}$として，レナード・ジョーンズポテンシャルの関数形は以下のようになります。
$$U(r_{ij}) = 4β \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^{6} \right] {\tag {13}}$$
$\beta$と$\sigma$は原子の種類に依存するパラメータです。
$r^{-12}$の項は核間距離が小さいときに生じる交換相互作用による強い反発力を，$r^{-6}$の項は双極子-双極子相互作用(ファンデルワールス相互作用)の主要項による引力を表していると考えることができます。
このとき，原子$j$が原子$i$に及ぼす力のベクトル$f_{ij}$は，
$$f_{ij}=\frac{48\beta}{r^2} \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\frac{1}{2}\left(\frac{\sigma}{r_{ij}}\right)^{6} \right]\bf r_{ij} {\tag {14}} $$
となります。
プロトタイプコードは以下のとおりです。
(ここでは周期境界条件と力の到達距離のカットオフを定めています。それに伴い，鏡像力の評価等をしています。これらについての説明は本稿で割愛しますが，このプログラムに関連した記事の投稿の際に説明する予定です)
"""
ビリアルの計算
Sept. 2017
"""
import numpy as np
import random
def sign(a,b): #周期境界条件のための設定
if b >= 0:
return abs(a)
else:
return -abs(a)
def virial():#ビリアルWを返す関数
#レナードジョーンズポテンシャルのパラメータ
beta=1
sigma=1
rcut = 3#カットオフ
r2cut = rcut**2
W=0 #ビリアル
fx[:] = 0
fy[:] = 0
fz[:] = 0
for i in range(N-1):
for j in range(i+1, N):
dx = x[i] - x[j]
dy = y[i] - y[j]
dz = z[i] - z[j]
# 近接 イメージ相互作用
if (abs(dx) > 0.5*L):
dx = dx - sign(L,dx)
if (abs(dy) > 0.5*L):
dy = dy-sign(L,dy)
if (abs(dz) > 0.5*L):
dz = dz-sign(L,dz)
r2 = dx**2+dy**2+dz**2
if (r2 < r2cut):
if (r2 == 0):
r2 = 0.0001
invr2 = 1/r2
wij = 48*beta*(sigma**12*invr2**3-0.5*sigma**6)*(invr2**3) # レナード・ジョーンズの原子間力
fijx = wij*invr2*dx
fijy = wij*invr2*dy
fijz = wij*invr2*dz
fx[i] += fijx
fy[i] += fijy
fz[i] += fijz
fx[j] -= fijx
fy[j] -= fijy
fz[i] -= fijz
W += wij
return W
# サンプル
N=100 #原子数
x = np.zeros([N],float)
y = np.zeros([N],float)
z = np.zeros([N],float)
fx = np.zeros([N],float) #力の配列の格納
fy = np.zeros([N],float)
fz = np.zeros([N],float)
L = int(N**(1/3)) # シミュレーションボックス: LxLxL
for i in range(N): # N原子をLxLxLのボックス内にランダムに配置する
x[i]=L*random.random()
y[i]=L*random.random()
z[i]=L*random.random()
virial()
ランダムな結果になります。5回ほど実行してみました。4桁から6桁までの値をとります。
231147.37139342312
6481.7477137863552
16588.246153314412
343999.10888889519
41240.556050538886
相互作用がない無い場合(理想気体)と比べると，ビリアルが正(負)だと相互作用により圧力は増加(減少)します。
参考文献
久保亮五 編, 大学演習 熱学・統計力学, 修訂版, 裳華房, (1998).