アルゴンのポテンシャルエネルギーの揺らぎから求めた臨界点

(1)

法政大学情報メディア教育研究センター研究報告 Vol.35 44 法政大学情報メディア教育研究センター研究報告 Vol.35 （2020）

原稿受付 2020 年 3 月 15 日

1．はじめに

気液臨界点は状態式が既知であれば，次の式のよ うに圧力 p の体積 V による 1 次微分と 2 次微分が 0 の点として定めることができる［1］． ∂ ∂    = ∂∂     = p V p V T T 2 2 0 （1）しかし状態式を計算から定めることは容易ではない．そこで通常の分子動力学シミュレーション（MD）［2］によって体積一定の条件のもとで，ポテ ンシャルエネルギー（Ep）の揺らぎの温度変化を求 める．温度が下がって，安定領域から不安定領域に入れば構造は大きく揺らぐのでその反映としてポテンシャルエネルギーが大きく揺らぐ．この変化から Spinodal 線［3］を決めることができる．Spinodal 線は安定領域と不安定領域の境界だからである［3］． Spinodal 線において最も高温の点を臨界点と定めることができる．この論文ではこの方法で臨界点を近似的に定める．分子動力学シミュレーションのアンサンブルは体積と温度一定の NTV である．この方法は一般的に安定領域と不安定領域の境界を定めるのに利用できると期待される．NTV-MD で は温度 T は指定された値の近傍で揺らぐが，今回は 毎ステップごとに速度をスケールしているため，温

アルゴンのポテンシャルエネルギーの揺らぎから求めた臨界点

Critical Point of Argon from Fluctuation of Potential Energy

片岡洋右1）2）

Yosuke Kataoka

1）_{法政大学情報メディア教育研究センター} 2）_{法政大学生命科学部環境応用化学科}

Critical point of argon was estimated from fluctuation of potential energy. Lennard-Jones potential was assumed in molecular dynamics simulation. Ensemble was canonical. Spinodal line was obtained from the temperature-dependence of the fluctuation of potential energy. Critical temperature was assigned as the maximum temperature in the spinodal line in the temperature-density space. The results were compared with the reported equation of state.

Keywords : Critical Point, Spinodal Line, Molecular Dynamics Simulation, Lennard-Jones Potential

度の揺らぎ幅は小さい．また体積 V は指定された値で一定である．このためこのアンサンブルは統計力学におけるカノニカルアンサンブルに相当する．このため指定された温度体積において構造の揺らぎが大きいかどうかをポテンシャルエネルギーの揺らぎから知ることができる．

2．分子動力学シミュレーション

アルゴンの臨界点と Spinodal 線を求めるため NTV 分子動力学シミュレーション（NTV-MD）を行った．使用したアルゴンの Lennard-Jones パラメータは表 1 に示した［4］．ポテンシャル関数 u(r) は次の式（2）に 示した．分子動力学シミュレーションの条件は表 2 に示した．式（2）におけるεはエネルギーパラメータであり，ポテンシャルの深さに相当する．またσはサイズパラメータであり原子直径に相当する． Lennard-Jones ポテンシャルで多くの液体−気体系を表すことができるので，具体的な計算はアルゴン表 1 アルゴンの Lennard-Jones パラメータ Table 1 Lennard-Jones Parameters for Argon［4］

ε/J σ/m 1.7258E-21 3.4282E-10

(2)

法政大学情報メディア教育研究センター研究報告 Vol.35 45 について実行するが，この論文の方法は多くの液体−気体系に適用できる． u r r r ( ) = _ _ − _ _        4ε σ 12 σ 6 （2）

3．Spinodal 線

密度が 0.5 g/cm3_{におけるポテンシャルエネルギー} の平均値 <Ep> とその標準偏差 MSD （Ep）の温度変 化を図 1 と図 2 に示した．この密度を調べたのは凝縮相の密度から van der Waals 式の場合の経験則を適用すると，近似的な臨界密度がこの値になるためである．物理学では揺らぎの大きさが問題になるがこ れは標準偏差の 2 乗である．またこの時の圧力 p の 温度変化を図 3 に示した． 図 2 からこの密度では，170 K 付近から MSD （Ep） が大きくなり始めるので，170 K 以上で安定状態であり，以下では不安定状態と読み取れる．こうした方法で求められる境界の温度には±5 K 程度の不確 定さがある．図 1 の <Ep> の温度変化と矛盾しない． 図 3 からは Spinodal 点付近では圧力の温度変化になんら際立った傾向の変化は現れないことが読み取れる．

4．原子配置

上の結論を視覚的に確認するため，図 4 と図 5 に原子配置の温度依存性を示した． これらの図から T＝135 K においては，構造の揺ら ぎが大きいことが分かる．一方 T＝170 K では構造の 揺らぎが 135 K の時ほどは大きくない．図 2 から 135 K におけるポテンシャルエネルギーの標準偏差が 170 K の時と比べ大きいことが確認できる．表 2 分子動力学シミュレーションの条件 Table 2 Conditions in Molecular Dynamics Simulation

item value

initial conﬁguration FCC number of molecules 864

ensemble NTV

velocity control Nose method boundary condition periodic

cut-oﬀ distance half of cell length

time step dt 1 fs

number of steps 1,000,000 potential function Lennard-Jones

application SCIGRESS-ME［4］

-6.5

-6

-5.5

-5

-4.5

-4

-3.5

-3

100

120

140

160

180

200 (Ep/1e-18)/J

(<

E

p

>

/

1 e

-18

)/

J

T/K

d=0.5g/cm

3 図 1 密度が 0.5 g/cm3_{におけるポテンシャルエネルギー} の平均値 <Ep)> の温度変化

Fig. 1 <Ep> vs. T Plot at d＝0.5 g/cm3

0

5

10

15

100

120

140

160

180

200

(MSD(Ep)/1e-20)J

(M

S

D

(E

p)

)/

1 e

-1

8 )/

J

T/K

d=0.5g/cm

3 図 2 密度が 0.5 g/cm3_{におけるポテンシャルエネルギー} の標準偏差 MSD （Ep）の温度変化

Fig. 2 MSD (Ep) vs. T Plot at d＝0.5 g/cm3

図 3 密度が 0.5 g/cm3_{における圧力 p の温度変化} Fig. 3 P vs. T Plot at d ＝0.5 g/cm3

(3)

法政大学情報メディア教育研究センター研究報告 Vol.35 46

5．臨界点

密度が 0.6 g/cm3_{以上の場合のポテンシャルエネル} ギーの標準偏差の温度変化を図 6 に示した．また，密度が 0.5 g/cm3_{以下の場合は図 7 に示した．} これらの図から求めた Spinodal 線を図 8 に示し， Kolafa と Nezbeda の状態式［5］と比較した．今回得られた Spinodal 線は状態式の結果［5］と良く対応している．臨界温度は Spinodal 線における最高温度である．このようにして得られた結果を表 3 で状態式の結果と比較した．一般に臨界温度を定めるのは容易ではない．臨界点近傍における圧力の等温線を精度よく求めるには分子数を十分大きくとる必要がある．ここでは 864 個という小さな系で状態式の結果に近い臨界温度を得たという点で注目に値する．

6．まとめ

NTV アンサンブルでの MD から得られたポテンシャルエネルギーの標準偏差の温度変化から，安定領域から不安定領域へ変わる点（Spinodal）を定めた．Spinodal 線の最大温度として臨界温度を決めた．この方法は相互作用の詳細には依存せずに，一般的に利用できる方法である．表 3 アルゴンの臨界定数 Table 3 Critical Constants of Argon

MD Kolafa & Nezbeda［5］ Lofti et al［6］

Tc/K 170 167 163

dc/（g/cm3） 0.5 0.51 0.5

pc/atm 66 44 51

図 4 密度が 0.5 g/cm3_{，T＝135 K における原子配置の} 時間変化，左と右では時間が異なる

Fig. 4 Atomic Conﬁguration at T＝135 K, d＝0.5 g/cm3_, Left and Right Figures Correspond to the Diﬀerent Time

図 5 密度が 0.5 g/cm3_{, T＝170 K における原子配置の} 時間変化，左と右では時間が異なる

Fig. 5 Atomic conﬁguration at T＝170 K, d＝0.5 g/cm3_, Left and Right Figures Correspond to the Diﬀerent Time

図 6 密度が 0.6 g/cm3 以上の場合のポテンシャルエネルギーの標準偏差の温度変化

Fig. 6 MSD (Ep) vs. T Plot d＞0.5 g/cm3

2

4

6

8

10

12

14

16

18

80 100 120 140 160 180

200

0.08

0.1

0.2

0.3

0.4

0.5 (M

S

D

(E

p)

/

1 e

-2

0)

/

J

T/K

図 7 密度が 0.5 g/cm3_{以下の場合のポテンシャル} エネルギーの標準偏差の温度変化

Fig. 7 MSD （Ep） vs. T Plot d＜0.6 g/cm3

図 8 赤は MD による Spinodal 線 Fig. 8 Red is Spinodal Line by MD

(4)

法政大学情報メディア教育研究センター研究報告 Vol.35 47

参考文献

［1］片岡洋右．熱力学と分子動力学シミュレーションによる粒子系の相転移．法政大学，2018，博士論文．

［2］ Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids. Clarendon Press, Oxford, 1989.

［3］片岡洋右．分子動力学法による液化過程と spinodal 線．法政大学計算科学研究センター研究報告． 2019, 34, p. 9-15. ［4］ SCIGRESS-ME ．富士通．http://www.fujitsu.com/jp/ solutions/business-technology/tc/sol/sgme/summary/, （参照 2020-06-22）．

［5］ Kolafa, J.; Nezbeda, I. The Lennard-Jones fluid: an accurate analytic and theoretically-based equation of state. Fluid Phase Equilibria. 1994, vol. 100, p. 1-34. ［6］ Lofti, A.; Vrabec, J.; Fischer, J. Vapour liquid equilibria

of the Lennard-Jones fluid from the NpT plus test parti-cle method. Molecular Physics. 1992, vol. 76, p. 1319-1333.