地域密着型エリア広告配信リクルートの地域密着型広告ドコイク?アドネットワーク Ads by ドコイク?[無料でホームページを作成] [通報・削除依頼]
[無料でホームページを作成] [通報・削除依頼]

分子動力学 理解への道

-----------------------------------------------------------------------------
■ 基礎知識

◇ セルサイズの条件
• MDセルのサイズは、境界条件の影響を受けないように大きく設定する。ポテンシャルのカットオフ距離の2倍は最低限必要。

◇ 初期状態の設定
• 原子の初期速度は、ボルツマン分布に従った乱数で与える。平衡状態においては、原子の速度はボルツマン分布に従うことが統計熱力学から分かっている。速やかに平衡状態に到達するためにボルツマン分布を使う。
• 初期速度として、いきなり大きな速度(温度)を与えることは好ましくない。高温のシミュレーションを行う場合でも、最初は低い温度で始め、徐々に温度を上げていく。
• 初期速度や初期配置が原因で系全体が並進運動や回転運動を起こさないように、重心位置や速度を調整することも必要である。
• 自然な状態になるためにはある程度の緩和時間が必要であり、分子動力学では、十分な緩和計算を行ってから物性値の算出を行う。

◇ 差分法の時間刻みの設定
• 差分法の時間刻みΔtは、精度を確保しつつ計算時間を短くするように適正に設定する必要がある。
• 一般に、原子の振動周期の1/1000〜1/200あたりを目安に設定される。
• 水素などの軽い原子を扱う場合は、時間刻みを小さく(Δt〜0.1fs)とる必要がある。水素が少量含まれる場合は、水素原子のみ時間刻みを小さくするマルチタイムステップ法により計算時間の削減が行われる。
• 固体系のシミュレーションではΔt=0.5〜5fs程度が一般的。

◇ 物性値のゆらぎの取り扱い
• 分子動力学の物性値はゆらぎが大きいため、算出には十分な時間・空間平均が必要。とくい応力はゆらぎが大きい。

◇ 系の制御とアンサンブル
• 熱統計力学により、熱力学的独立変数は3つであるから、3つの状態変数を制御することができる。
• とくに、何の制御もしない場合は、粒子数N、エネルギーE、体積Vが保存される(一定に制御される)ミクロカノニカルアンサンブル(NEV)となる。下記に(NEV)と他のアンサンブルも含めて特徴を記す。
 ミクロカノニカルアンサンブル(NEV):エネルギーと体積が一定となる。温度と圧力のゆらぎは小さい。
 カノニカルアンサンブル(NTV):温度Tと体積Vが一定となる。NEVと比較して、圧力と全エネルギーとポテンシャルエネルギーのゆらぎは大きい。
 等圧アンサンブル(NPH):圧力Pを一定に制御する
 等圧・等温アンサンブル(NTP):温度と圧力を一定に制御する。圧力、体積が設定値のまわりを大きくゆらぐ。NEVと比較して、全エネルギーとポテンシャルエネルギーのゆらぎは大きい。
 等応力アンサンブル(NσH):応力σを一定に制御する
 等応力・等温アンサンブル(NσT):応力と温度を一定に制御する
※ Hはエンタルピー

◇ 原子間ポテンシャル
• 精度を最も上げるには、電子状態計算をベースとした第一原理分子動力学が理想である。しかし、計算系は100個程度に制限され、長い時間ステップの計算も難しいのが現状である。そのため、材料の特性に合わせてさまざまな原子間ポテンシャルが提案されてきた。

□ 2体ポテンシャル:計算時間が短いことや、ポテンシャル形状が単純で現象の解釈が容易であるなどの理由で用いられる。空孔や表面などの欠陥まわりのエネルギーは再現できないとされている。
 代表的な2体ポテンシャル:レナードジョーンズポテンシャル、モースポテンシャル、ジョンソンポテンシャル。

□ 3体ポテンシャル:欠陥(点欠陥、転位、粒界など)の解析で用いられる

• EAM系のポテンシャル:金属系で用いられる。EAM系として、EAM、MEAM、FS、GEAMが挙げられる。
 EAM:レナードジョーンズポテンシャルの2.3倍の計算時間がかかる
 MEAM:EAMの関数形状の改良、スクリーニング関数、角度項の導入などがされている。レナードジョーンズポテンシャルの20倍程度の計算時間がかかる
 FS:EAMポテンシャルと比べて関数形状が単純であり、コーディングしやすいことから多く用いられてきた。
 GEAM:2元系へ柔軟に対応できるポテンシャルとして広く使われている
 RGL:EAMの一種

• クラスターポテンシャル:角度依存の項が入ったものとしてSWポテンシャルがある。結晶構造が変化しない計算に向いている。
 SW:シリコンに適用される

• ボンドオーダーポテンシャル:低配位数(2〜4)の挙動を再現する。結晶構造が変化する計算に向いている。
 ターソフポテンシャル:シリコンと炭素に適用される
 ブレーナーポテンシャルもある。

• イオンポテンシャル:イオン結合を再現する。表面や界面など構造が急激に変化する場所では電荷が不均一に分散するため、電荷移動を考えた可変電荷ポテンシャルが必要となる。可変電荷ポテンシャルでは、MDステップごとに電荷の分布の計算が必要なため、通常の分子動力学のおよそ100倍以上の計算時間がかかるとされる。
 EAM+可変電荷ポテンシャル:NiOなどの金属系の酸化物に適応される
 ターソフ+可変電荷ポテンシャル:SiO2などの共有結合系の酸化物
 イオンポテンシャルを使う際には、周期境界条件の周期性の特徴を活かして、フーリエ変換により、遠方のコピーセルの影響を効率的に計算するエワルド法が用いられる。ただし、エワルド法は計算量が大きく、通常の分子動力学の10倍程度になると言われている。

• 分子内・分子間ポテンシャル:有名な分指導力学の汎用ソフトウエアのAMBERなどで提供されている。
 高分子などは、個々の原子間結合の挙動よりも分子鎖としての挙動に着目する場合が多い。基本的に最初に定義した分子内の結合は切れないと考える。したがって、全体のエネルギーは、分子内相互作用(分子内の化学結合のエネルギー)と分子間相互作用(分子間のファンデルワールス相互作用やクーロン相互作用)の和になる。

• 反応場(原子間ポテンシャル):結合の生成・解離を取り扱える。ReaxFFがその一例。
 ReaxFF:レナードジョーンズポテンシャルの250倍程度の計算時間がかかると推定される

◇ 原子間ポテンシャルについてのその他の事項
• 残念ながら、全ての物性値を再現する原子間ポテンシャルの開発は不可能であることが分かっている(第一原理分子動力学を使えばその問題は解決するが計算可能な原子数が通常数百のレベルに落ちてしまう。悩みどころではある)。
• 候補となるポテンシャルが選択できたら、それらの物性の再現性の調査を行う。
• 重要なことは、計算の目的を明確にして、計算の目的となる現象をポテンシャルが再現しているかという点を調査する。
• 目的の現象のみを再現できていればよいというわけではなく、凝集エネルギーや格子定数、弾性定数などのごく基本的な物性の再現も必要である。さもないと構造が不安定になったり、合金などの2元系のポテンシャルへの拡張がしにくくなったりする。
• 実績のあるポテンシャルのほとんどでは基本的な物性は再現されている。

-----------------------------------------------------------------------------
■ 分子動力学(MD)法とモンテカルロ(MC)法の使い分け
MD法: 支配方程式は古典力学の運動方程式、決定論
 流体のような原意位置の自由度が大きく、ダイナミックスや非平衡状態が重要な場合に用いる
MC法: 支配方程式は統計力学のボルツマン分布、確率論
 固体などのような原子位置が有限の自由度を持ち、途中のダイナミックより最終的な平衡状態が重要な場合に用いる
-----------------------------------------------------------------------------
■ 物性の評価に関する基礎知識

◇ 平衡原子間距離と線膨張係数
• 平衡原子間距離は、系の応力がゼロになるときの原子間距離である。
• 絶対零度では、ポテンシャル関数により解析的に求めることもできるが、有限温度においては、原子振動により原子間距離が常に変動するため、解析的には求められない。
• 分子動力学においては、平衡原子間距離を求めるために、まず、NTPアンサンブルにより、ゼロ応力となるセルサイズの平均値を求める必要がある。
• 温度が高くなればなるほど熱膨張するので、平衡原子間距離は大きくなる。
• 温度条件を変えたNTPアンサンブルによりセルサイズの温度依存性を計算することで、線膨張係数(線膨張率)を求めることができる。
• 線膨張係数には温度依存性があるため、興味のある温度の近くで計算する必要がある。

◇ 比熱
• 定積モル比熱は、温度Tに関する内部エネルギー(全エネルギー)Eの変化の割合である。
• 分子動力学で比熱を求める場合は、さまざまな温度において全エネルギーを求めて、温度に対する勾配を計算する。
• 固体の比熱は、デューロン−プティの法則により、元素によらず、3*気体定数となることが知られている。

◇ 輸送係数
• 輸送係数とは、物質拡散や熱の流速とその勾配を対応付ける係数であり、流速 = 輸送係数 * 勾配 という関係がある。拡散係数、熱伝導係数、粘性係数などがその例である。

◇ 動径分布関数
• 完全結晶の動径分布関数は、第1ピークから第2と第3と、近接原子に対応してピークが明瞭に観察される。
• アモルファス構造は、第1ピークは変わらないが、第2と第3ピークが統合され、中距離の秩序がなくなる。
• 上記のようなアモルファル構造の特徴は、最近接原子との距離は原子間距離によって定まるため結晶と変わらず秩序性があるが、構造の無秩序性のため、距離が離れるにつれ徐々に秩序を失っていくことに起因している。

◇ 反応
• 多くの物理現象は熱活性化過程であり、あるエネルギー障壁を熱エネルギーのアシストで乗り越えることにより現象が進んでいく。
• 拡散などの現象が起こる頻度(速度)はR(=6D/δ2)は、下記のように表すことができる。
 R = ν * exp(−ΔE/kBT)
 ここで、νは頻度因子、kBはボルツマン定数、Dは拡散係数、δは空孔のジャンプする距離である。
• エネルギー障壁が高い場合、障壁を乗り越える確率は低く、数千、数万振動周期の時間が必要となる。
• エネルギー障壁が高いと、拡散前の状態まわりで単純に熱振動している時間が非常に長く、まれにエネルギー障壁を乗り超えることになる。
 このような現象を再現するのに、分子動力学は効率的ではない。なぜなら、拡散が起こる前の熱振動状態が計算時間のほとんどを占めてしまうからである。エネルギー障壁が低ければ(<0.5eV程度)、十分な計算時間をかければ計算可能であるが、エネルギー障壁が高くなると現実的な計算時間では現象は起こらなくなり、分子動力学の適用が事実上不可能になってしまう。
 この解決方法としては下記の3つがある。

1. 温度を上げて現象を加速する。
 この方法は、融点以上に温度を上げることができず、温度を上げるとまったく異なる現象(相転移など)が起こる可能性があり、効果は限定的になってしまう。

2. 反応経路解析
 ある原子系の状態を原子座標の3N次元の状態空間上の一つの点と考え、反応(拡散)前の状態の点と反応(拡散)後の状態の点の間にいくつかの中間状態(レプリカと呼ばれる)を設定して、仮想的な弾性バネで結合する。弾性バネで結合することにより、各状態の点は一定の距離を保ちながら反応経路上に配置されることになる。ただし、初期に設定した配置が最小エネルギー反応経路を表しているわけではないので、反応経路と垂直なバネ力を最小化するようにバネを少しずつずらすことによって、最もエネルギーが低い経路を探索する。このような手法は、ナッチド・エラスティックバンド法(NEB法)と呼ばれ、材料系の反応経路解析においては標準的に用いられている。
 実際の計算では、レプリカの設定が大きな課題となり、現象に合わせてさまざまな手法が提案されている。このレプリカの設定手法とNEB法(およびその改良)を合わせて反応経路解析と呼ぶ。
 反応経路解析は、分子動力学の時間スケール(〜ns)では起こらない現象を推測する有効なツールである。

3. 加速分子動力学法
 反応経路解析では活性化エネルギーΔEは求められるが、頻度因子νを求めることはできないため、実際にどの程度の頻度で拡散などの現象が起こるかはわからない。すなわち、現象の起こる頻度Rを求めることはできない。
 一方、分子動力学の計算では、活性化エネルギーが低い場合に限られるが、現象が実際に起こるのでRを求めることができる。また、Rの温度依存性のデータの傾きから、活性化エネルギーを求めることができる。しかし、ΔEが高い場合は、分子動力学では現象がまったく起こらないという問題が生じる。そこで、分子動力学を加速するためのさまざまな手法が提案されてきた。フォルターが提案したハイパーダイナミックス法は、原子間ポテンシャル関数Vに何らかのブーストポテンシャルΔV(s)を追加して、現象を速やかに起こすことを目的としている。ここで、sは何らかの原子構造の幾何学的パラメーターである。
 Vb = V + ΔV(s)
 ΔEが高い状態では、深いエネルギー曲面の底に系がトラップされている。そこで、エネルギー曲面の底の部分にブーストポテンシャルΔV(s)を追加して、系の状態を深い底から追い出すことを考える。このブーストポテンシャルは任意に設定可能であり、通常は、何らかの原子構造の幾何学的パラメータsを用いることが多い。たとえば、結合長さが平衡結合長を保たなくするようなバイアスをかけることで、拡散などの原子結合を切る現象を加速することができる。これをボンドブースト法と呼ぶ。
 ただし、ブーストポテンシャルが乗り越える先の状態への経路の鞍点(サドルポイント)にかかってしまうと、起こる現象の反応経路と活性化エネルギーを変えてしまうことになるので、その範囲に配慮が必要となる。
 ブーストポテンシャルが設定できたら、分子動力学計算を行う。ポテンシャルがかさ上げされている分だけ、早く現象が起こる。これは本来ならば、ポテンシャル井戸の底の部分に滞在する時間をスキップしていることに相当する。計算上は、現在の状態にかさ上げされているΔVを使って、1MDステップの分子動力学の時間刻みΔtMDが、Δtboostに相当すると解釈する。時間を強制的に進める作業をMDステップごとに行う。
 Δtboost = ΔtMD * exp(ΔV/kBT)
 加速分子動力学は低温でも計算可能である。高温では分子動力学と同じ結果が得られている。加速分子動力学ではハイパーダイナミックス法のほか、メタダイナミクス法などのいくつかの手法が提案され、それぞれ発展している。

□ 温度・圧力制御の注意点
 MC法ではアンサンブルに応じて基準となるエネルギー表式が決まるため問題はないが、MD法は運動方程式を追跡するだけなので本質的に外界と熱のやり取りのないエネルギー保存系の計算であり、いわばNEVアンサンブルしか扱えない。実際に計算結果を実験データと比較するためには温度と圧力をコントロールする必要があり、そのための処方が提唱されている。ただ、それらの処方は必然的に外界との熱のやり取りが必要となるため、運動方程式とは別の操作もしくは新たな支配方程式が加わることになり、その点がMD法の計算結果の信頼性を損ねる一因となっている。
 系の温度は原子の運動エネルギーの総和から常に計算できる。そこで系の温度を設定温度に変えたい時には原子の速度を調整すればよい。その処方としては、
(i) 全原子の速度を一斉に定数倍する(速度スケーリング法 [TP1])
(ii) 温度を上げたい時には加速、下げたい時には減速する粘性係数的な仮想変数を導入し、その仮想変数の支配方程式は設定温度と系の温度との差に比例して時間変化数する(Nose-Hoover法 [TP1, 2])
(iii) 原子の速度を設定温度のボルツマン分布のものとランダムにおきかる
などがある。どれも一見乱暴なやり方であるが、系が熱平衡に達した時にはどの処方でも(ある誤差の範囲内で)正しい熱平衡分布を与えることが示されている。[TP4]
 圧力制御に関しても似た状況で、周期境界条件を課し、そのセル内での圧力は熱力学のビリアル定理から常に計算できるため、設定圧力と異なる場合には、それに合わせて基本セルの大きさを変えればよい。処方として、
(i) セルの大きさを微少量ずつ変化させる(セルスケーリング法)
(ii) セルのサイズや形そのものを仮想変数として、セル自身の運動方程式は設定圧力と系の圧力との差に比例して形状が変化する形を取る(Anderson法 [TP3]、Parrinello-Rahman法 [TP5])
などがある。温度・圧力制御ともに仮想変数の変化に伴うオーバーシュートや非現実な相変化が起こりうるため、MD法において非現実的な結果が出やすい部分であるため注意が必要である。
[TP1] L. V. Woodcock et al., Chem. Phys. Lett., 10 (1971) 257.
[TP2] W. G. Hoover et al., Phys. Rev. A, 31 (1985) 1695.
[TP3] H. C. Anderson et al., J. hem. Phys., 72 (1980) 2384.
[TP4] S. Nose et al., J. Chem. Phys., 81 (1984) 511.
[TP5] M. Parrinello and A. Rahman, Phys. Rev. Lett., 45 (1980) 1196.
-----------------------------------------------------------------------------
■ プログラミングで重要になる知識

◇ ポテンシャルのカットオフ
• カットオフは、計算結果に影響が出ないように十分に遠くで滑らかに切り捨てる必要がある。単純に切り捨てると、カットオフ距離においてエネルギーの不連続性が生じてしまうため、エネルギー値がカットオフ距離でゼロとなるように関数をシフトする必要がある。これをシフトポテンシャルと呼ぶ。
• 単にシフトするだけでは、エネルギー値は連続であっても、エネルギーの微分値(力)が連続ではなくなる。そのため、エネルギーの微分値がカットオフでゼロになるフォースシフトポテンシャルがよく用いられる。
• ポテンシャル関数にカットオフ関数をかけることによって滑らかに切り捨てる方法もよく用いられる。

◇ ブックキーピング法
• 分子動力学では、ポテンシャルのカットオフ距離内の相互作用の計算のみを行う。このため、どの原子がカットオフ距離内にあるかの台帳(ブック)が計算ステップごとに必要となる。
• カットオフ距離より少し大きい範囲に存在する原子の番号を台帳に記憶しておき、台帳の外の原子がカットオフ距離内に入ってくる可能性がない間は台帳の更新を行わないことにより、計算量を大幅に削減できる。このような手法をブックキーピングと呼ぶ。
• ブックキーピング法は分子動力学の高速化に欠かせない手法である。しかし、台帳の更新間隔、台帳の範囲などを適正に設定する必要がある。

◇ 差分法
• 分子動力学では、計算精度を高めるためにさまざまな差分法が提案されている。
• 速度ベルレ法:精度が比較的高く計算負荷が低いことからよく用いられる差分法
• ギア法:精度を上げるために、加速度より高次の時間微分量などを用いて、次のステップの状態を決める。加速度以上の高次の微分量の扱いは容易でないため、過去の時刻の加速度を用いる手法が一般的である。ギア法は、ギアの予測子ー修正子法とも呼ばれ、次のステップの予測値をもとに、運動方程式より修正量を計算し、次のステップの修正値を定めている。ギア法は計算負荷は高いが、時間刻みが小さい場合は、ベルレ法より精度が高いことが分かっている。
-----------------------------------------------------------------------------
■ Alの多結晶
 数10nm程度の粒子サイズを境にホール・ペッチ則と逆ホール・ペッチ則が入れ替わることや、亀裂先端から発生した転位が結晶粒界に到達した後、粒界が新たな転位発生源となることを見出している。
[1] Shimokawa et al.,
[2] 日本金属学会関東支部『ここまでできる。計算材料科学入門』
[3] ホール・ペッチ則 (http://tkkoba.blog.fc2.com/blog-entry-80.html, http://kakiken1984.blog103.fc2.com/blog-entry-146.html)
[4] 逆ホール・ペッチ則 (http://www.nims.go.jp/apfim/nanostrength.html)
-----------------------------------------------------------------------------
参考文献
[1]泉聡志ら『実践分子動力学シミュレーション』森北出版
[2] 日本金属学会関東支部『ここまでできる。計算材料科学入門』
-----------------------------------------------------------------------------
QRコード
携帯用QRコード
アクセス数
ページビュー数