カー・パリネロ法
出典: Wikipedio
|
|
カー・パリネロ法（カーパリネロほう、Car-Parrinello method、略してCP法とも言う）は、1985年、カー（R. Car）とパリネロ（M. Parrinello）によって考案された手法である[1]。従来のバンド計算で用いていた、行列要素の対角化（固有値問題）を用いずに固有値（及び固有ベクトル）を求めることにより、計算を大幅に高速化することを実現した。これにより、系の電子状態と共に、その構造の最適化（この部分は古典的分子動力学法を用いる）も可能とした。
目次 |
名称に関して
本手法は、第一原理分子動力学法(First-Principles Molecular Dynamics:FPMD)、またはab initio MD(MD: Molecular Dynamics)などとも呼ばれる。それぞれ微妙に異なるものを意味している場合もある。現在では、カー、パリネロによるオリジナルな手法が用いられることはほとんどなく、より効率化、高速化を図った手法に置き換わっている。ただし、これら現在主流となっている手法も広い意味でのカー・パリネロ法の範疇にあると言える。
カー・パリネロ法使用の利点
- 計算量の減少 : 対角化による手法（固有値問題）を使用しないので、計算量のオーダーをN3から最大 NlogN 程度まで減らすことができる。N は使用する基底関数の展開数。
- メモリの節約 : 同様にして必要なメモリも大体N2から、NMのオーダーにすることが出来る。M はバンドの数。ここで基底関数の数は、バンド数より一桁以上大きいことが前提（←平面波基底の場合、N >> M）
- 系の構造の分子動力学計算や最適化が電子状態計算と同時に行える。
ラグランジアンが出発点
<math> L = \sum_{i,\vec{k}} \int m |\dot \psi_{i,\vec{k}} (\vec{r}) |^2 d\vec{r} + \sum_I {1 \over 2} M_I \dot\vec {R_I}^2 - E_{tot}[{\psi_{i,\vec{k}} },{\vec R_I }] + \sum_{i,j,\vec{k}} \Lambda_{i,j,\vec{k}} (\int {\psi^*}_{i,\vec{k}} (\vec{r}) \psi_{j,\vec{k}} (\vec{r}) d\vec{r} - \delta_{i,j}) </math>
- <math>i</math> : バンドの指標
- <math>\vec{k}</math> : k点
- <math>m</math> : 波動関数に対する仮想質量
- <math>\psi</math> : 波動関数
- <math>M_I</math> : イオン芯（原子核）部分の質量（Iは指標）
- <math>\vec{R}_I</math> : イオン芯の座標の位置ベクトル
- <math>E_{tot}[]</math> : 系の全エネルギー
- <math>\Lambda_{i,j,\vec{k}}</math> : ラグランジュの未定係数
ドットは時間微分を表し、上付きの"*"は複素共役を表す。
<math> \int {\psi^*}_{i,\vec{k}} (\vec{r}) \psi_{j,\vec{k}} (\vec{r}) d\vec{r} = \delta_{i,j} </math> 規格直交性（基底関数に対する制約）
<math> {d \over dt} ({ {\delta L} \over {\delta \dot{\psi^*}_{i,\vec{k}} } }) = {{\delta L} \over {\delta {\psi^*}_{i,\vec{k}} } }</math> 波動関数部分
<math> {d \over dt} ({ {\partial L} \over {\partial \dot{\vec{R}_I} } }) = { {\partial L} \over {\partial \vec{R}_I } } </math> ← イオン芯部分
以上から、以下の二つの運動方程式が得られる。tは時間（波動関数に対する時間tは仮想的な時間で実時間ではない）、δは変分（汎関数微分）を意味する。Hはハミルトニアン。
<math> m \ddot{\psi}_{i,\vec{k}} (\vec{r},t) = - { {\delta E_{tot}} \over {\delta {\psi^*}_{i,\vec{k}} (\vec{r},t) } } + \sum_j \Lambda_{i,j,\vec{k}} \psi_{j,\vec{k}} (\vec{r},t) = - H \psi_{i,\vec{k}} (\vec{r},t) + \sum_j \Lambda_{i,j,\vec{k}} \psi_{j,\vec{k}} (\vec{r},t) </math>
<math> M_I \ddot{\vec{R}_I} = - \nabla_{\vec{R}_I} E </math>
参照 : ラグランジュ力学
上記、最後の2式から系の電子状態及び系の構造の最適化を行う。波動関数に関しての仮想的な運動方程式は、時間の2階微分を1階微分に置き換えると最急降下法となる。他にもいくつかの手法（例 : 共役勾配法など）への発展形がある。
過去、現在、未来
カー・パリネロ法が出てくる以前は、電子状態の計算（バンド計算）を行いつつ、構造の最適化（〔準〕安定構造の探索）を同時に行うことは、計算（対角化手法←固有値問題）で求められる原子間に働く力から、手動でユニットセル内の原子を動かして次のステップに回すか（手動で動かす段階で一旦計算が終了してしまう）、計算量が膨大な対角化手法をそのまま用いて、計算から得られる力をもとにユニットセル内の原子を分子動力学を用いて動かすことが行われていた。いずれにしても、大変効率が悪く、扱える原子数はせいぜい数個のオーダーであった(*)。
カー・パリネロ法の出現は、この扱える原子数を一気に数十個のオーダーに引き上げた。数十個という規模ならスラブ近似を用いれば比較的扱い易いシリコンの表面系の安定構造を求める計算や、構造最適化の動力学的な過程を追っていくことが可能となった。
(*)そもそもカー・パリネロ法が登場する以前は、通常のバンド計算（電子状態計算）で構造の最適化（〔準〕安定構造の探索）を同時に遂行することはほとんどなく、表面などのような複雑な系での構造は、実験によって観測された値を流用することがほとんどであった。
初期のカー・パリネロ法では、系の原子（より正確にはイオン芯または原子核）を分子動力学によって解くのと同時に、電子の波動関数を仮想的に時間発展していくものとして、波動関数に関する仮想的な運動方程式を分子動力学的手法と同じ形式を用いて原子系と連動して解いていた（上記、ラグランジアンの式変形参照）。これは後により効率の良い方法へと発展していくこととなる。
当初は、専ら擬ポテンシャル＋平面波基底による方法で、カー・パリネロ的手法が実現されたが、その後、全電子手法であるAPW法に対し、カー・パリネロ的手法を取り入れたものが出現している[2,3]。更に、タイトバインディング法との結合 (TBMD) や、混合基底を使った手法でもカー・パリネロ的手法を導入したものが登場している。
電子状態を解くための手法面でも、初めは電子の波動関数を最急降下法やベレの方法で逐次的に解く方法が使われたが、その後、共役勾配法（いくつかの改良版がある）や、より洗練された方法（例 : RMM-DIIS法、ダビソン法）が使われるようになっている。
波動関数を計算する部分（電子状態計算部分）も、原子の動力学と同時に解くのではなく、原子を動かす時間幅を大きくとるようにして、原子を動かす毎に電子状態部分を常にボルン‐オッペンハイマー面まで収束させるようにする手法が主流となっていった（ペインのアルゴリズム[4]）。また、この手法では原子の移動（構造の最適化）も同時に行われるが、この時原子の移動量による波動関数の変化を外挿により予想する方法[5]など、高速化、効率化を図る数多くの提案がなされている。
最近では、オーダーN法や、ハイブリッド法内での利用などの拡張も行われている。
温度に関して
ユニットセル内の原子（イオン芯）を分子動力学手法で扱う場合は有限温度での計算が可能であるが、一方、同時に行う電子状態計算部分は密度汎関数法の枠内での計算なので絶対零度での計算となる。ペインのアルゴリズムのように、原子を動かす毎に電子状態をボルン‐オッペンハイマー面に収束させる場合、計算上の支障はないが、電子状態計算はあくまで絶対零度での計算であることに注意が必要である。電子状態計算部分も有限温度へ拡張するアプローチも存在する。
カー・パリネロ法で必要な近似、手法、道具
- 上述のように、最急降下法、共役勾配法、RMM-DIIS法などを使って、波動関数の更新（電子状態部分の計算）を行う。
- 断熱近似 : 電子状態の計算と共に、ユニットセル内の原子の分子動力学計算（→原子を動かす）を行うので、この近似が成立しない系には適用できない。
- ヘルマン-ファインマン力 : 原子（イオン芯）の分子動力学計算を行うためには、原子間に働く力を電子状態計算から求める必要がある。
- 圧力（ストレス） : 初期の頃は、ユニットセル内の原子（イオン芯）の構造の最適化のみが行われたが、後に圧力やストレスも計算してユニットセルの大きさや形そのものも最適化の対象となった。つまり電子状態の計算と同時に、ユニットセルの内部構造及びユニットセル自身の最適化（〔準〕安定構造の探索）も現在では行われるようになっている。
- グラムシュミットの直交化 : 電子状態計算において、基底は直交していなければならない（←少なくともセルフコンシステントな計算が収束した段階では）。このため直交化（手法←グラムシュミットの方法以外の直交化手段もある）が必要。
- スーパーコンピュータ、超並列マシン : カー・パリネロ法の出現により、計算速度、効率は飛躍的に上がったが、電子状態部分の計算には依然として大量の計算資源が要求され、2003年においても原子数が100個を越えるような大きな系の計算の実現のためには、スーパーコンピュータや超並列マシンが必須である。
- ワークステーション、パーソナルコンピュータ(PC) : カー・パリネロ法が登場した当初はPC上での第一原理分子動力学計算などは夢のような話であったが、同法が利用されるようになってから20年近くが経った2003年では、手法そのものの向上及び、ワークステーションやPC等の飛躍的な性能向上も相まって、スーパーコンピュータのような巨大な計算資源を使わなくとも中規模程度（〜数十原子からなる系）の第一原理分子動力学計算なら（←計算条件、目的にも依る）、ワークステーションや数台のPCクラスターあるいは単独で計算遂行が可能になっている。
参考文献
- [1] R. Car and M. Parrinello, Phys. Rev. Lett., 55 (1985) 2471.
- [2] J. M. Soler and A. R. Williams, Phys. Rev. B40 (1989) 1560.
- [3] J. M. Soler and A. R. Williams, Phys. Rev. B42 (1990) 9728.
- [4] M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, Review of Modern Physics, Vol. 64, No. 4, (1992) 1045. （総合的な解説記事としても重要）
- [5] T. A. Arias, M. C. Payne and J. D. Joannopoulos, Phys. Rev. B45 (1992) 1538.
