どぼるざーくとゆく、高分子シミュレーション

lammpsやOCTAなどの高分子シミュレーターについて四苦八苦しながら学んでいくブログです。Twitterアカウント@lammps_octa

【お勉強】分子動力学法-1(原子に作用する力)

はじめに(2020/01/12加筆)

今回の記事はlammpsで主に使われる分子動力学法の解説をしたいと思います。

(今回の記事は学問的なものなので問題点があればご指摘頂けると幸いです。)

 

 

 

簡単な説明

分子動力学法(MD,Molecular Dynamics)とは全ての原子の運動について運動方程式を解くことにより、その軌跡を追跡する計算方法です。運動方程式とは『ニュートンの第二法則』の運動方程式です。

f:id:toal_peg_simu:20200106095111p:plain ニュートンの運動方程式 - Wikipedia

 

 この方程式にシミュレーション開始時の『初期位置』と『初期速度』、各時間における『原子に作用する力』を代入すると、次のstepにおける原子の位置と速度を求めることが出来ます。

 

原子に作用する力 is 何?

先ほどの運動方程式に影響する『原子に作用する力』には以下のようなものがあります。

シュレーディンガー方程式(⇒第一原理分子動力学)

f:id:toal_peg_simu:20200106101524p:plain  シュレーディンガー方程式 - Wikipedia

シュレディンガー方程式を解くと電子の波動関数が分かり、そこから『フェルマン・ファイマン力』を求めることが出来ます。この力は原子核に作用するため、原子の運動に影響を与えます。

電子状態が変化するような系(レーザーによる励起など)を取り扱う際に有効です。ただし時間に依存するシュレディンガー方程式は解くことが難しいので、いくつかの拘束条件が必要になります。この拘束条件を系に合わせて設定しなければシミュレーション結果に悪影響が生じます。

 

②ポテンシャル関数(⇒古典分子動力学)

原子に作用する力を『ポテンシャル関数』で近似してシミュレーションをします。有名なポテンシャル関数には『Lennard-Jonesポテンシャル』があります。

U(r)=4\epsilon \left[\left({\frac  {\sigma }{r}}\right)^{{p}}-\left({\frac  {\sigma }{r}}\right)^{{q}}\right] レナード-ジョーンズ・ポテンシャル - Wikipedia

 うまくポテンシャルを設定すればシュレディンガー方程式よりも計算量を減らすことが出来るため、大きな系についても扱う事ができるようになります。一方でポテンシャルを事例に応じて調べる必要があります。

前回の記事でシミュレーションでは古典分子動力学を使用しました。

 

原子の位置と速度を決定する(運動方程式の数値積分)

原子の位置と速度を決定するためには運動方程式に当てはめて解く必要があります。パソコン上で微分方程式を解く場合には差分方程式に変換します。

この差分方程式を解くにはいくつかの手法があります。

Leap frog法

{\displaystyle {\begin{aligned}x_{i+1}&=x_{i}+v_{i}\Delta t+{\frac {1}{2}}a_{i}\Delta t^{2},\\v_{i+1}&=v_{i}+{\frac {1}{2}}(a_{i}+a_{i+1})\Delta t.\end{aligned}}} リープ・フロッグ法 - Wikipedia

二次近似で微分方程式を解く手法である。この手法のメリットは時間可逆性があるということです。また運動のエネルギーが保存されるという性質ももちます。

ただしこの手法では加速度が速度に非依存であり、時間ステップを変えることができないので注意が必要です。

 他にはverlet法などがあります。(解説したかったですが、イマイチ違いが分からないのでまた今度にします)

 

Leap frog法(詳しく)

 Leap frog法の所で『時間可逆性』『エネルギー保存則』のところがイマイチ分からなかった、というコメントをもらったので少し調べてみました。

 

時間可逆性

まず時間可逆性が守られる意味についてなんですが、そもそも現実世界は非・時間可逆性らしいです。分かりやすい資料があったので以下に載せておきます。

https://www.gakushuin.ac.jp/~881791/materials/Irreversiblity09.pdf

https://en.wikipedia.org/wiki/Time_reversibility

ここで注目すべきことは『ニュートンの力学法則で成立している系では、時間の向きが可逆である』ということです。つまり分子動力学法ではニュートンの第二法則が系を支配しているため、時間の可逆性が成り立つことになります。

ここでもう一度数値積分について考えてみましょう。シミュレーションにルンゲクッタ法を用いると、そもそもルンゲクッタ法に時間可逆性がないため、シミュレーション自体も(ニュートン力学の系なのに)時間可逆性が失われてしまいます。このような論理の破綻をなくすために、分子動力学法では時間可逆性がある『Leap frog法』などが使われているのかもしれませんね。

 

エネルギーの保存

修正前の記事では時間可逆性とエネルギーの保存に因果関係があるかのように書いていましたが、実際は独立の性質でした・・・(すみません、間違っていたので書き直しました)

恐らくここで言いたかったのは、各粒子に初期値として与えられた運動エネルギーが長時動シミュレーションをしても変動しないことがleap frog法のメリット、ということが言いたいのではないでしょうか。

 

 

 

さいごに

いかがでしたか。

今回はlammpsでよく使う分子動力学法の簡単な説明と運動の推測手法について説明を行いました。

難しい表現で分からない部分があればコメントなどで指摘してください。