Lennard-Jones流体の分子シミュレーション(2)分子動力学vsメトロポリス法 -

以前、分子動力学の実装を作ったけど、諸条件、諸パラメータをこれと合わせることで、メトロポリス法と分子動力学で結果を比較できる。メトロポリス法を使ったことがなかったので、やってみようと思った。

分子動力学では、NVE一定の計算が、一番基本的になるのに対して、メトロポリス法では、NVT一定の計算が、一番基本的になる。メトロポリス法でも、NPT一定の計算をやる方法は考えられているし、分子動力学でも、barostatやthermostatを使えば、NVTやNPT一定の計算ができることにはなってるけど。また、メトロポリス法に不利な点として、輸送係数を計算する方法がないと思う。

全部実装して計算してから、気付いたけど、以下に書いたようなことは、Verletが、1967年の論文でやっていた。この論文は、Lennard-Jones流体の分子動力学計算を報告した最初の論文だそうで、この論文にさらっと書かれた数値積分法は、その後Verlet法と呼ばれるようになった

Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules
https://doi.org/10.1103/PhysRev.159.98

メトロポリス法の論文は1953年に出ている。この論文は、与えられた確率分布に従う乱数列を生成するメトロポリス・ヘイスティングス法の最初の例が記述された論文として、物理以外の分野にとっても、重要な論文となった。

Equation of State Calculations by Fast Computing Machines
https://doi.org/10.1063/1.1699114

何故か、Wikipediaにまでページがある
Equation of State Calculations by Fast Computing Machines
https://en.wikipedia.org/wiki/Equation_of_State_Calculations_by_Fast_Computing_Machines

著者は5人いるけど、その一人Rosenbruthが書いた

Genesis of the Monte Carlo Algorithm forStatistical Mechanics
https://doi.org/10.1063/1.1632112

によると、統計力学の問題をターゲットにしようと決めて、中心となるアルゴリズムを考えたのは、Rosenbruthっぽい。また、最初は、普通に"分子動力学"をやろうと思ったが、計算機のスペックが厳しそうなので別の方法を考えたともある。論文で使われた計算機は、MANIACのようで、メトロポリスが開発を指揮したらしい。

メトロポリスらの論文で実際に計算してるのは二次元剛体球模型だけど、殆どの部分は、模型の詳細とは関係ない。ついでなので、論文を眺めた。

古典統計力学では、アンサンブル平均の計算が基本的になる。アンサンブル平均は、相空間での積分
$\bar{F} = \left( \displaystyle \int F(\mathbf{q} , \mathbf{p}) \exp(-E/k_{B} T)d^{M}p d^{M}q \right) / \left( \displaystyle \int \exp(-E/k_{B} T)d^{M}p d^{M}q \right)$
のこと。 $k_{B}$ はボルツマン定数、 $E$ はエネルギー(ハミルトニアン)で、相空間上の関数。 $F$ は、相空間上の任意の関数。 $T$ は温度。相空間の次元は $2M$ で、空間次元が3、粒子数がNなら、 $M=3N$ であるが、メトロポリスは、空間次元が2で考えているので、 $M=2N$ 。例えば、 $F=E$ の場合は、アンサンブル平均は、内部エネルギーを与える。

よく出てくるアンサンブル平均の計算では、運動量pに関する積分は解析的に実行して消去できることが多いので、位置に関する積分のみが残る。例えば、ちょっと計算すれば分かるけど、内部エネルギーUは、
$U = \dfrac{M k_{B} T}{2} + \left( \displaystyle \int e^{-\phi / (k_{B} T)}d^M q \right) / \left( \displaystyle \int \exp(-E/k_{B} T)d^{M}p d^{M}q \right)$
のように書ける。相互作用項がない場合は、理想気体の結果として、統計力学の教科書に書いてあると思う。

運動量が動く範囲は力学的に決まるが、位置座標の範囲は、人間が決める。特に理由が無い限り、立方体領域にするのが簡単なので、普通、そうすると思う(周期境界条件を課すことが多いので、正確にはトーラス)。立方体の体積Vと、粒子数N,温度Tを外部パラメータとして指定することになるので、正しくカノニカル・アンサンブルの計算を行うことになる。

計算したいものは単純だけど、粒子数が100程度でも、M=200とかM=300なわけで、低次元で有効な数値積分のアルゴリズム(台形公式やガウス求積)を使って評価するのは難しい。しかし、ポテンシャルエネルギーが大きい配置(例えば、二粒子が極めて近接してるような)は、積分に小さな寄与しかしないので、現代で言う棄却サンプリングによって、粒子の配置を $\exp(-\phi /k_{B} T)$ に比例する確率で発生させ、このような配置に対して、 $F(\mathbf{q})$ の平均を取れば、アンサンブル平均を効率的に計算できる。

ところで、適当な統計力学の教科書を開くと、hをPlanck定数として、分配関数
$Z = \dfrac{1}{N!h^{M}} \displaystyle \int \exp(-E/k_{B} T)d^{M}p d^{M}q$
を計算すれば、自由エネルギーとか、エントロピーとか圧力とかの熱力学諸量は、分配関数の対数を取ったり微分したりする数学的操作で全部計算できるよ〜と書いてある(今の場合、分配関数Zは、N,V,Tの関数)。が、メトロポリス法では、分配関数そのものは計算できないので、こういう話は全然役に立たない。統計力学の教科書は破り捨てよう。

メトロポリス法で計算できるのは、アンサンブル平均なので、熱力学量をアンサンブル平均に帰着する必要がある。メトロポリスらの論文では、圧力を扱っている(論文のタイトルにもある通り、状態方程式を調べたいというのが動機だったので)。圧力は、統計力学の教科書的な定義では
$P = \dfrac{1}{\beta Z} \dfrac{\partial Z}{\partial V}$
である。この式を形式的に計算すると、ハミルトニアンを体積で微分した関数のアンサンブル平均という、よく分からないものが得られる。体積は、Hamiltonianの変数ではないので、素朴な意味では、この式は意味をなさない。

(メトロポリスらの論文では使われてないけど)たまに見かける議論として、体積の代わりに、Hamiltonianの位置と運動量のスケール変換(正準変換になるように、位置をr倍したら運動量は1/r倍する)を考えて、r=1での微係数を計算をすると、ビリアル圧力が得られる。正確には(空間次元を3とすれば) $V=r^3 V_{0}$ に対して、 $\dfrac{\partial}{\partial V} = \dfrac{1}{3r^2 V_{0}} \dfrac{\partial}{\partial r}$ なので、体積と空間次元で割る必要がある。形式的な議論に見えるけど、配位空間内の積分領域でなく、相空間内の積分領域に作用する一径数(正準)変換群を考えて、この変換群で、分配関数がどう変化するのか見ていて、やってることは正当である。正準変換で、相空間の測度が不変に保たれることが重要。

ビリアル圧力は、位置のみの関数で書けるので、メトロポリス法で計算するのに適している。分子動力学でも、ビリアル圧力は、瞬時圧力を計算するのに一般的に使われている。分子動力学の場合は、マクロな圧力は、瞬時圧力の時間平均として得る。

スケール変換に限らず、もっと一般に座標の線形変換を、相空間の正準変換に拡張できる。具体的には、 $A \in GL(n,\mathbf{R})$ に対して、 $\left(\begin{matrix} A & 0 \\ 0 & (A^T)^{-1} \end{matrix} \right) \in Sp(2n,\mathbf{R})$ を全粒子に一様に施す。 $A$ が単位行列の定数倍であれば、さっきのスケール変換の状況。 $GL(n,\mathbf{R})$ には、直交変換群 $O(n)$ が含まれているが、これは、系全体を回転する(あるいは鏡に映す)だけなので、なんの影響も及ぼさない。直感的には、この線形変換は、粒子のいる立方体領域を変形する操作に相当する。

無限小線形変換で考えると、
$r_{\alpha} \mapsto r_{\alpha} + \sum_{\beta} \epsilon_{\alpha \beta} r_{\beta}$
の形で、歪みテンソル $\epsilon_{\alpha \beta}$ は対称としていい(反対称成分は、無限小回転に相当)。Hamiltonianにこの変換を施して、歪みテンソルで微分すると、応力テンソル(virial stress)の成分が計算できる。virial stressは、運動量流束密度と見ることもでき、粘性係数を計算するのにも使われる。

他に得られる物理量としては、内部エネルギーの熱的揺らぎから、定積比熱が得られる(はず)。内部エネルギーやエントロピーは、直接測定できる量ではないが、定積比熱は測定できる量だけど、多原子分子だと、分子内振動の自由度からの寄与が正しく評価されない気がする(これは古典統計力学では、どうやっても計算できないので)。単原子分子でも、分子間相互作用からの寄与が、どれくらい正しく評価されるものか分からない。

エネルギーの二乗のアンサンブル平均は、 $\beta=1/(k_{B}T)$ として、以下のように計算できる
$\bar{E^2} = \dfrac{M^2+2M}{4\beta^2} + \left(\dfrac{M}{\beta} \displaystyle \int \phi e^{-\beta \phi} dq + \displaystyle \int \phi^2 e^{-\beta \phi} dq \right) / \left( \displaystyle \int \exp(-\phi/k_{B} T) dq\right)$

エントロピーは、相空間上の関数のアンサンブル平均として直接計算する方法はないと思う。エントロピーが得られると、ヘルムホルツの自由エネルギー $F=U-TS$ がわかり、 $Z=e^{-F/k_{B}T}$ によって、分配関数が計算できてしまう。

以上は模型によらない一般論だけど、2次元剛体球(rigid sphereと書いてるけど、現在は、hard sphereの方が、通りがいいと思う。2次元なので、circleあるいはdiskだけど、rigid circleとかhard diskとは言わない)の計算では、もう少し、固有の工夫が要る。

剛体球模型は、粒子同士が重なると、ポテンシャルエネルギーが無限大、重ならなければ、ポテンシャルエネルギーは0というモデルなので、 $\exp(-\phi /k_{B} T)$ は、0か1しか取らない。従って、棄却サンプリングでは、重なる配置は即棄却、そうでない配置は即採用と比較的アルゴリズムが単純になる。

一方、ビリアル圧力の計算には、粒子に働く力を知る必要があり、これは、通常の場合、ポテンシャルエネルギーの微分で計算されるが、剛体球模型では、許可される配置に於いて、ポテンシャルエネルギーは常に0であって、そのまま適用できない。剛体球模型で力が働くのは、粒子同士が衝突した一瞬であって、メトロポリスの論文では、この衝突力の評価を、動径分布関数(ある粒子から一定距離にある粒子密度)から評価する式を与えている。動径分布関数も、解析的に決めることはできないが、数値的に外挿することで評価したと書いてある。

あと、メトロポリス法では、乱数生成が必要だけど、疑似乱数生成アルゴリズムに何を使ったかは書かれてない。当時は、von Neumann考案のmiddle-square methodと、Lehmerの方法というものが知られていたっぽい。Lehmer乱数は、線形合同法の特殊ケースでもある。多分、どっちかの方法を使ったんだろう。現在は、どっちも使われてないアルゴリズム。

二次元剛体球模型は、実験との比較ができないtoy modelなだけでなく、計算機が速くなった今となっては、本質的じゃない部分で却って面倒が多い。というわけで、普通に、3次元のソフトポテンシャル模型を考え、ポテンシャルのモデルとしては、ポピュラーなLennard-Jonesポテンシャルを使ってメトロポリス法を実装した。希ガス気体は理想気体と変わらなすぎるので、希ガス液体(今回は希ガスアルゴン)の計算をやる。で、冒頭に書いたように、分子動力学や実測結果と比較する。

分子動力学の計算では、全粒子対で相互作用計算をすると、粒子数Nに対して、O(N^2)なので、カットオフ距離を設定して、計算量をO(N)にしている。メトロポリス法は、最初に一回だけ全粒子対の計算をするのみなので、カットオフ距離は別に必要ない。ただ、今回は、条件を揃えて計算するために、カットオフ距離を設定できるようにした。

Lennard-Jonesポテンシャルで計算するのに、カットオフを用いる場合、内部エネルギーや圧力に、カットオフで削った部分の長距離補正項を足す必要がある。 $g(r)$ を動径分布関数とすると、内部エネルギーと圧力は
$U = \dfrac{3}{2}Nk_{B}T + 2\pi N \rho \displaystyle \int_{0}^{\infty} \phi(r) g(r) r^2 dr$
$p = \rho k_{B} T - \dfrac{2 \pi \rho^2}{3} \displaystyle \int_{0}^{\infty} \phi'(r) g(r) r^3 dr$
などと書ける。 $\rho$ は粒子数密度。長距離補正は、カットオフ距離より遠い部分の積分で
$U_{LRC} = 2 \pi N \rho \displaystyle \int_{r_c}^{\infty} \phi(r) g(r) r^2 dr$
$p_{LRC} = -\dfrac{2 \pi \rho^2}{3} \displaystyle \int_{r_c}^{\infty} \phi'(r) g(r) r^3 dr$
という形。Lennard-Jones流体に対する補正項の計算では、 $g(r) \approx 1$ と近似するのが一般的。こうすると、積分が計算できる。圧力の長距離補正は、密度の二乗で効いてくるので、気体の時はともかく、液体では、かなり大きな影響を持ちうる。

ついでに、一般のvirial stressに対する長距離補正は
$\sigma^{ij}_{LRC} = \dfrac{\rho^2}{2} \displaystyle \int_{|\mathbf{r}| \gt r_c} \frac{r_i r_j}{|\mathbf{r}|} \phi'(|\mathbf{r}|) g(|\mathbf{r}|) d\mathbf{r}$
という形。
$p_{LRC} = -\dfrac{1}{3}(\sigma^{xx}_{LRC} + \sigma^{yy}_{LRC} + \sigma^{zz}_{LRC})$
が成り立つ。。

カットオフを入れると、ポテンシャルエネルギーが不連続になるので、微分すると、力が発散するのを気にする人もいるらしい。エネルギーが連続になるように修正したポテンシャルエネルギーを使うこともある。
cf)Numerical Experiments on the Stochastic Behavior of a Lennard-Jones Gas System
https://doi.org/10.1103/PhysRevA.8.1504

Verletは気にしなかったようなので、これは、気にしなくてもいい気がする(というか、入れてない)

アルゴンのLennard-Jonesパラメータを使って、粒子数10万、密度1407(kg / m^3)、温度85(K)あたりで、液体と思われる平衡状態を生成して、これを初期状態に、カットオフ距離を変えつつ、分子動力学とメトロポリス法で圧力を計算した結果及び圧力補正項は以下のようになった。rcutがカットオフ距離で、sigmaはLennard-Jonesパラメータのσ。通常、rcutはsigmaの2.5〜10倍程度の範囲で設定するっぽい

rcut / sigma	p(atm)@MD	p_LRC(atm)	T(K)@MD	p(atm)@MC
2.5	1.93	-312	85.26	-2.57
5.0	9.80	-39.1	85.02	11.98
7.5	11.25	-11.6	84.95	12.55
15.0		-1.45		12.90
30.0		-0.181		10.43

メトロポリス法の温度は、85(K)で、積分点を100万個サンプリングした。カットオフ距離が大きくなると、カットオフによる計算高速化の意味がなくなるので、分子動力学では、大きなカットオフはやってない。圧力の計算値は、長距離補正込みで、p_LRCは圧力の長距離補正。今は液体を扱ってることもあって、標準的な範囲のカットオフ距離だと、長距離補正の影響は、結構でかい。

Verletの論文では、rcut/sigmaは、2.5や3.3を使ったと書いてある。Verletの論文を読んだのが、この計算した後だったせいで、Lennard-Jonesポテンシャルのパラメータは、Verletの使った数値とはかすかに異なる(例えば、私はσを3.41Åにしたけど、Verletは、3.405Å)

Verletの論文で、rcut/sigmaが小さめなのは、粒子数が少ないせいもある。私が、粒子数10万を選んだのは、特に根拠があったわけでもないけど、今回くらいの条件で、長距離補正が支配的にならないようにしたいと思うと、粒子数10万で、多すぎるとも言えないように思う。

で、実験と比較したい。以下は、やや古い論文ではあるが、測定結果に基づいて、液体アルゴンの温度と圧力から密度を計算するフィッティング式が書かれている。

Density of liquid nitrogen and argon as a function of pressure and temperature
https://www.sciencedirect.com/science/article/abs/pii/0031891460900422

85(K)での測定結果はないが、この式に従って計算すると、
85(K),1(atm)での密度は1409(kg/m^3)
85(K),10(atm)での密度は1411(kg/m^3)
85.25(K),1(atm)での密度は1407(kg/m^3)
85.25(K),10(atm)での密度は1410(kg/m^3)
85.5(K),10(atm)での密度は1408(kg/m^3)
86(K),10(atm)での密度は1405(kg/m^3)
などとなるので、計算結果は、全くの的外れってことはなさそう。

最初は、(温度は85Kで)、1400(kg/m^3)、1405(kg/m^3)などで計算したが、いずれも負の圧力が実現されたので、次に1407(kg/m^3)を選んで計算をしたら、正の圧力が得られたので、それから文献調査をした。同じ温度で、1(atm)と10(atm)でも、密度が2~3(kg/m^3)くらいしか変わらない割に、そこそこの値を出していると思う。

Lennard-Jonesポテンシャル自体が、そもそも雑な近似(パラメータの決め方も正しい方法とかないし)なので、実験結果との一致がしょぼくても仕方ないけれど、もう少し、ちゃんと測定結果との比較をしてみようと思って、以下の論文の結果と対応する計算結果を出した。論文では、温度と密度と圧力が与えられてるので、温度と密度を指定して、圧力を計算する。

Measurement and correlation of the (pressure, density, temperature) relation of argon I. The homogeneous gas and liquid regions in the temperature range from 90 K to 340 K at pressures up to 12 MPa
https://doi.org/10.1006/jcht.1994.1048

T(K)	density(kg/m^3)	p(atm)@exp	p(atm)@MD	T(K)@MD	p(atm)@MC
90	1381.34	9.65	25.97	90.20	26.81
100	1320.44	19.58	25.59	100.05	22.50
110	1250.17	19.76	10.49	109.68	10.42

粒子数は10万。
MDのカットオフ距離は、σの5倍。
メトロポリス法のカットオフ距離は、σの30倍で、積分点は100万点。

分子動力学で初期平衡状態を作って、NVE一定の分子動力学とメトロポリス法による計算を行った。とりあえず、圧力を高い精度で予測するには至ってない。

今回の計算に使用したコードを、gistに置いておく。
https://gist.github.com/vertexoperator/49737c13e432bf13e26c851683a6821d
パラメータとかは、mainの中でハードコードしてるのを手動で書き直して、使ってたので、使い勝手はよくない。

差し当たって、適当に検索して出てきた数値を使ってシミュレーションした。LJパラメータは、文献によって、そこそこバラツキがあるので、仮に、たまたま良いパラメータを引いて、ある範囲の条件下で、熱力学量がよく計算できても、パラメータいじれば合うよねという感じがしなくもない。なので、Lennard-Jonesパラメータにも根拠が欲しい。結局、実験的に決めるより、第一原理的に計算する方法があればベスト思われ、それには量子力学が必要だけど、随分昔にやった

量子力学はLennard-Jonesポテンシャルを再現できるか
https://m-a-o.hatenablog.com/entry/20121217/p2

を眺めると、いまいちLJパラメータの計算値と実験値は合わないと書いている。今適当に検索すると、このにっきを書いた翌年に

Theoretical Investigation of the Dispersion Interaction inArgon Dimer and Trimer
https://doi.org/10.1016/j.procs.2013.05.247

という論文が出てて、CCSD(T)/aug-cc-pVTZで計算した場合、Ar-Arの平衡分子間距離が3.75Åあたりになると書いてある。普通のLJポテンシャルで、エネルギー最小点が一致するようにパラメータを決めると、σが3.34Åくらいに相当。この論文で引用されてた

A highly accurate interatomic potential for argon
https://doi.org/10.1063/1.466051

を見ても、平衡分子間距離は、ほぼ同じ値。Verletの使ったσが、3.405Å。εについては、Fig1(a)のグラフを見ると、平衡分子間距離におけるエネルギーが、(無限遠に置いた時を基準として)約-100(/cm)=-0.0124(eV)。上のシミュレーションで使った値は、eV換算で、0.01034(eV)。

LJポテンシャルの式の形は単純すぎるので、別の関数形を使うことも考えるべきなのかもしれない。安直にパラメータを増やすのはよくないが、双極子の次は、四重極も...と考えるのは自然なことで、実際に"精密化学"界隈では、r^(-8)やr^(-10)に比例する項も計算されるのが普通らしい。適当に検索した所、London分散力の延長として、四重極を考慮した人は、1938年には存在している
Quadrupole Contributions to London's Dispersion Forces
https://doi.org/10.1063/1.1750184
古典分子動力学で、そのようなポテンシャルを使っていけない理由は別にない。

このへんの問題は、おいおい追求するとして、とりあえず、CCSD(T)/aug-cc-pVTZ計算の追試のため、Gamessに以下のような入力ファイルを食わせた所、核間距離は、3.8008895(ANGS.)で、TOTAL ENERGYは、-1054.1027355965(hartree)となった。論文より平衡分子間距離が少し大きいけど、まぁ許容範囲と思う

! rungms Ar2.inp  > Ar2.log
 $contrl runtyp=optimize CCTYP=CCSD(T) NUMGRD=.T. ISPHER=1 coord=unique $end
 $system TIMLIM=20000 MEMORY=10000000 $end
 $basis gbasis=ACCT $end
 $guess guess=Huckel $end
 $data
argon dimer
C1
 Ar  18 0.0 0.0 0.0
 Ar  18 0.0 0.0 3.7
 $end

同じようにして、希ガス二量体に対するCCSD(T)/aug-cc-pVTZレベルの計算をした結果を、まとめると、以下のようになった。

	Rm(Angs.)	Em(hartree)	E0(hartee)	(E0-Em)/k_B(K)	R_nn(Angs.)	融点(K)	沸点(K)	融解熱/R(K)
Ne	3.0916867	-257.6254676649	-257.6252957319	54.293	3.13	24.56	27.07	40.29
Ar	3.8008895	-1054.1027355965	-1054.1023022117	136.853	3.76	83.80	87.30	141.92
Kr	4.1146791	-5504.4997943600	-5504.4992507992	171.644	4.01	115.79	119.93	197.25

Rm:平衡分子間距離(計算値)
Em:平衡分子間距離におけるエネルギー(計算値)
E0:核間距離を100Åにした時のエネルギー(計算値)
k_B:ボルツマン定数
R_nn:希ガス固体の最隣接原子間距離(Kittelの教科書に載ってるらしいのを孫引き)
融点:出典Wikipedia(多分、1気圧の数値？)
沸点:出典Wikipedia(多分、1気圧の数値？)
融解熱:出典Wikipedia
R:気体定数