magattacaのブログ

日付以外誤報

雰囲気量子化学入門(後編) ~ハートリー・フォック法の課題と密度汎関数理論〜

雰囲気量子化学入門・後編!前回に引き続き、明日知ったかぶれる量子化学を目指して用語をひろっていきます!

前編ではシュレーディンガー方程式からハートリー・フォック方程式までを辿りました。 後編ではハートリー・フォック法の問題点とその解決策、そしてDFTについてお勉強!

f:id:magattaca:20210123183315p:plain

全体の流れ

まず、全体の流れを再掲します。だいたいこんな感じの話です。

f:id:magattaca:20210123012656p:plain

今回は下半分、ハートリー・フォック方程式の誤差から出発したいと思います。

ハートリー・フォック法の誤差(電子相関)

さて、シュレーディンガー方程式から近似で導いたハートリー・フォック方程式で得られるエネルギーは、全体のエネルギーの99%以上を見積もることができます *1。残る誤差は1%程度ですが、この大きさは化学に適応する点で非常に問題となります。

なぜかというと、化学で考慮したい課題は、化学反応の活性化エネルギーやイオン化エネルギーといった問題ですが、これらは全エネルギーの絶対値ではなく、 そのエネルギー差(相対値)によって結果が左右されます。このエネルギー差の大きさが、ハートリー・フォック法の誤差と同程度となるため、 化学的精度にとって無視できない誤差となってしまいます。*2

f:id:magattaca:20210123172606p:plain

この誤差は電子相関と呼ばれ、レフディンにより電子相関エネルギーは以下のように定義されています。

f:id:magattaca:20210123012924p:plain

大まかにいうと、電子相関は「多電子系において1個の電子の位置を定めると、他の1個の電子の存在確率が第1の電子の位置の影響を受ける」 というようなことを意味するようです。*3

ではなぜハートリー・フォック法ではこの相関を一部しか取り込めず、誤差となってしまうのでしょうか?要因と解決策を辿ってみましょう。

電子相関の起源とポスト・ハートリー・フォック法

ここではより具体的な電子相関の起源と、その解決策として考案されたハートリー・フォック法の延長線上にある手法を辿ります。

ざっとこんな感じです。

f:id:magattaca:20210123172706p:plain

動的電子相関と静的電子相関

まず、電子相関は大きく動的電子相関静的電子相関の2種類に分類されます。*4

動的電子相関はハートリー・フォック法が独立粒子モデルに基づいていることに由来します。 このモデルでは「1つの電子に着目した平均場近似」となっているため、電子間の位置の相互作用平均的なものに置き換えています。 つまり、2つ以上の電子の衝突散乱といった効果が含まれません。この、「平均化により表現できなくなった部分」が動的電子相関の起源となっています。

格好いい表現でいうと、「ハートリー・フォック波動関数相関カスプ条件を満足しないため、動的電子相関が欠如しており、 短距離電子反発を余計に取り込んでしまう」(???)ということになるようです。

要するに、実際には電子同士は互いに反発するのでゼロ距離まで近づく確率はとても低いけど、ハートリー・フォックはそんな状況が存在する可能性も含めてしまう*5。余分な反発まで考えてしまうので、結合エネルギーの過小評価や、結合長の過大評価につながったりしてしまう、ってことでしょうか??

f:id:magattaca:20210123172856p:plain

一方、静的電子相関は1つのスレーター行列式(電子配置)を使って波動関数を近似していることに由来します。

格好いい表現では「いくつかの電子状態(電子配置)がエネルギー的に近く、擬縮退している場合にみられる 配置間の長距離相互作用(擬縮退効果)」(???)となります。

良くわからないのでもう少し具体的に、例として水素分子を水素原子に解離することを考えます。

この時、閉殻ハートリー・フォック(RHF)法では「1つの空間軌道に2つの電子をペアで詰める」ので、2つの水素原子の分子軌道のうち、一方に電子2個、もう一方に電子が0個のイオン化した状態しか得られません。ですが、実際にはそれぞれ1個ずつ電子を持つ水素原子として解離した状態の方がエネルギー的に安定です。従って"不自然"な状態を過大評価してしまうことになっています。

つまり、「1つよりも同時に2つを考えた方がいい結果が得られる時があるのに、それが上手くできないので問題だ!」って感じでしょうか?

f:id:magattaca:20210123173022p:plain

・・・良くわかりませんがとにかく誤差があるんだ!!!どげんかせんといかん!

ポスト・ハートリー・フォック法

ハートリー・フォック法で取りこぼしてしまう電子相関を考慮するために、様々な手法が考案されています。

上図に示したように大きく次の3つの分類が可能です。

  • ① 動的電子相関を考慮するためのポスト・ハートリー・フォック法
  • ② 静的電子相関を考慮するための多配置 (multi-configuration) SCF法
  • ③ どちらも考慮する多参照 (multi-reference)電子相関法

このうち、ポスト・ハートリー・フォック法では、仮想軌道を用いることで励起状態の電子配置(励起配置)を考慮します。 原理的には「全ての励起配置を考慮すれば、用いている基底関数の空間内で、電子相関を厳密に取り扱える」ことになります。 一方で、励起配置の数が増えるほど計算コストが増大するという課題があります。

f:id:magattaca:20210123013306p:plain

ポスト・ハートリー・フォック法には大きく① 配置間相互作用(CI)法、② 摂動法、③ クラスター展開(CC)法があります。 それぞれこんな感じです。

f:id:magattaca:20210123013432p:plain

CI法は線型結合で表現されており変分法による近似でCI係数を決定します。

一方、摂動法はその名の通り近似に摂動法をもちいています。 ハミルトニアン無摂動ハミルトニアン摂動項分割し、できるだけ少ない逐次近似で真の解に近づくような分割の仕方が目指されます。

これらに対して、CC法では演算子指数関数となっている点に特徴があります。CI法と同じ数の変数でより多くの励起配置を考慮でき、 より高い精度の結果が得られると言う利点があります。

いずれも「電子の励起配置をどこまで考慮するか?」によって呼び方が変わります。

MP2とかCCSD(T)とかなんとなく聞いたことありましたけどポスト・ハートリー・フォック法だったんですねー。

多配置SCF法、多参照電子相関法

ポスト・ハートリー・フォック法以外も少し触れておきます。

まず、静的電子相関を考慮する方法として、多配置SCF法があります。 基本的考えは「複数のスレーター行列を使って波動関数を表現する」というもので、配置関数によるCI展開を用います。 配置関数の選び方に複数あり、代表的なものはCASSCF法です。*6

また、両方(動的+静的)の電子相関を考慮する方法として多参照電子相関法があります。 多配置SCF法で静的な部分を考慮しつつ、ポスト・ハートリー・フォック法のいずれかを用いて動的な部分を取り込む、組合せのアプローチです。*7

以上、ハートリー・フォック法の拡張を見てきました。これらはいずれも「波動関数に基づいて電子相関を考慮するアプローチ」ですが、計算コストが大きくなるという問題点があります。

これに対して、より低い計算コストで電子相関を取り込むことができる手法があります。それが密度汎関数理論(DEnsity Functional Theory, DFT)です。

ようやく本題!

密度汎関数理論

DFTでは、波動関数を求めるのではなく、代わりに「電子密度から出発して電子相関を考慮」します。 この「波動関数電子密度」というのがDFTの最大の特徴です。

これまで、解けない「波動関数に関するシュレーディンガー方程式」を解くため、近似を頑張ってハートリー・フォック方程式へと導きました。 対して、DFTではシュレーディンガー方程式と等価な「電子密度に関するコーン・シャム方程式」を用います。

どうしてこれがOKになるのか、順番に見ていきましょう!

お話の概要を再掲しときます。

f:id:magattaca:20210123013602p:plain

ホーヘンベルグ・コーンの定理とエネルギー汎関数

まず、ホーヘンベルグ・コーン(Hohenberg-Kohn)の定理と呼ばれる基本定理から出発します。

以下の2つの定理からなります。

f:id:magattaca:20210123013621p:plain

① 存在定理は要するに、「系のすべての物理量基底状態の電子密度を与えれば一意的に決まる、すなわち電子密度の汎関数である」ということを意味するそうです。

上図に記載の通り、エネルギーも電子密度の汎関数(エネルギー汎関数)となっています。右辺はそれぞれ、T[ρ]:運動エネルギーV_ext[ρ]:外部ポテンシャルエネルギーV_ee[ρ]:電子間相互作用エネルギーの項を表します。

また、② 変分原理は、「種々の電子密度 (ρ)のうち、エネルギー汎関数(E[ρ])に代入して得られるエネルギーが最小となるような電子密度(ρ_0)が、 系の基底状態の電子密度であり、「正解」である」といった意味になるそうです。つまり、「電子密度(ρ)を色々と変化させて、とにかく条件に合うものを見つけてしまえば解けたことになる」、と! *8

たぶん凄いことなんですが、私の能力ではイマイチ凄さがつかみきれない!悔しい!

要するに、「等価だから、シュレ-ディンガー方程式を解く代わりにDFTのアプローチを採用した方が計算楽かもよ!」ってことですね。知らんけど

エネルギー汎関数の形

さて、等価な問題と書きましたが、正確な電子状態が求まるのは、「エネルギーと電子密度の間の対応関係を正確に表す汎関数がわかっている場合」に限られます。

しかし、残念ながらそのような汎関数の形はわかっておらず、結局DFTにおいても近似に基づいた汎関数を使うことになります。

近似として以下のようなものが提案されましたが、精度が低く*9、実用性に問題がありました。

f:id:magattaca:20210123013723p:plain

これらに対し、コーンとシャムが導入した近似により、DFTは実用上有用な方法となりました。

コーン・シャム近似

コーンとシャムが導入した近似はだいたい以下のようなものです。

f:id:magattaca:20210123013932p:plain

厳密な表式がわからず近似の精度の悪さに関係していた運動エネルギーに着目し、正確に表せる部分とそれ以外の残りの部分に分けました。 前者はコーン・シャム軌道と呼ばれる軌道を導入して表し、残りの表しきれない部分を交換相関汎関数に押し込めてしまいました。

また、加えて、古典的な電子間クーロンエネルギーで近似できない量子力学な電子間相互作用についても、この交換相関汎関数に押し込められました。

f:id:magattaca:20210123014025p:plain

「運動エネルギーの”正確に表せる部分”」にもう少し踏み込みます。ここでは、現実の複雑な電子同士の相互作用のある系の代わりに、「相互作用のない仮想的な系(Kohn-Shamの補助系」を導入します。この補助系は「現実の系と同じ電子密度(+同じ全エネルギー)」をもつように設定されており、 補助系の運動エネルギー汎関数は厳密に表現することができます。 *10

こうして得られたコーン・シャム法のエネルギー汎関数では、独立粒子運動エネルギー項(T_KS)と 原子核-電子間ポテンシャルエネルギー汎関数(V_ne)、古典的な電子-電子間のクーロン相互作用によるエネルギー項(J)、 そして交換相関汎関数(E_XC)とに分割されて表されるようになりました。

交換相関汎関数の形

さて、コーン・シャム近似でよくわからない部分を押し込められた交換相関汎関数ですが、やはり厳密な形がわからず、 様々な近似で表すことがアプローチされています。

こんな感じ。

f:id:magattaca:20210123014113p:plain

いよいよB3LYP出てきました!!

さて、交換相関汎関数の最も単純な近似は一様電子ガスモデルから導出された「局所密度近似(local density approximation、LDA)」でした。 この精度を高めるために、電子密度の勾配(∇ρ)の効果(電子密度の変化に対するエネルギーの依存性)を考慮してLDAを補正した、「一般化された密度勾配近似(generalized gradient approximation, GGA)」が考案されました。

しかし、依然として「交換エネルギーをポテンシャル関数で表す近似」では化学的精度を得るには粗い場合があり、密度汎関数法とハートリーフォック法とのハイブリッドと呼ばれる方法が提案されました。この「混成GGA汎関数(hybrid-GGA)」では、ハートリー・フォック交換項を使ってGGAに補正が加えられています。*11

交換相関汎関数の具体的な形として多くが提案されていますが、通例「交換汎関数+相関汎関数」と両者を組み合わせた名称で呼ばれます。 これは、「電子のスピンに由来する効果を記述する交換汎関数」と「電子相関に由来する効果を記述する相関汎関数」とにわけて それぞれ開発されてきた経緯によります。

で、Beckeによる交換汎関数(B)とLee-Yang-Parrによる相関汎関数(LYP)の組み合わせに、ハートリー・フォック交換項を加えた混成GG汎関数B3LYPであったと!なるほど!!

コーン・シャム方程式とその解き方

さて、交換相関汎関数についてもわかりコーン・シャム法のあらましがわかってきました。あとはこの問題をどう解けば良いか?、です。

が、ちょっと用語が多すぎて目標を見失なってしまったので、もう一度振り返ります。

とにかく我々の目標は「系(分子)の状態(物理量)」を知ることです。シュレーディンガー方程式では波動関数を求めればOKでした。 対して、DFTでは基底状態の電子密度を求めます。で、この電子状態を求めるのに必要な(だけど正確な形がわからない)汎関数の近似のために、 コーン・シャム軌道や交換相関汎関数が導入されました。ここで、コーン・シャム軌道は相互作用のない独立電子系(コーン・シャムの補助系)を考えているので。 1つの電子に関する軌道です。

1電子軌道とエネルギーの問題・・・前回の記事でみたハートリー・フォック法と同じ形式の問題ですね!

ハートリー・フォック法における解法のアプローチはこんな感じでした。

f:id:magattaca:20210123014223p:plain

コーン・シャム法においても同じアプローチを用いることができます。

変分原理よりエネルギーを最小化するものが最も良い近似なので、コーン・シャム軌道の規格直交条件を付加条件としてラグランジュの未定乗数法を使うことで、コーン・シャム方程式と呼ばれる方程式が導出されます。

こんな感じ。

f:id:magattaca:20210123014257p:plain

導いたコーン・シャム方程式は、ハートリー・フォック方程式と同じ形の1電子方程式となっているので、同様に解くことができます。

ハートリー・フォック・ローターン方程式を導いた時のように、基底関数系を用いた展開によって行列問題に変換(コーン・シャム行列)すれば、 SCFの手続きによる繰り返し計算によって解くことができます。

前回の記事で「行列問題への置き換えはDFT計算においても主流の解き方」と書きましたが、確かにコーン・シャム法も同じアプローチで解けましたね!

また、この置き換えにおいて、DFTでも基底関数系が用いられました。前回見たように基底関数系の例として有名なものにPople型基底関数(ex. 6-31G)があるのでした。

これで論文で良く見る表現B3LYP/6-31G」が「交換相関汎関数/基底関数系」の組み合わせを表すことがわかりましたね!!これで今日から知ったかぶれる!!

コーン・シャム法の振り返りと解の意味するもの

無事、密度汎関数理論において解くべき方程式(コーン・シャム方程式)とその解法まで辿ることができました。最後に、以上の手続きで得られたものについて振り返りたいと思います。

まず、密度汎関数理論の導入としてハートリー・フォック法におけるエネルギーの誤差の問題(電子相関)から出発しました。 ハートリー・フォック法の延長にあるポスト・ハートリー・フォック法では、電子相関を「基底状態励起状態との電子配置間の相互作用」として取り込みました。

対して上図で示したDFTのコーン・シャム方程式において、交換相関ポテンシャルは「交換相関エネルギー(E_XC)を電子密度(ρ)で微分したもの」として定義されています。つまり、コーン・シャム法においては、電子相関はポテンシャル汎関数として取り込まれていることになります。

ハートリー・フォック法とコーン・シャムDFTの式を見比べてみます。

f:id:magattaca:20210123014421p:plain

ほとんど同じ形ですが、コーン・シャムDFTではハートリー・フォック法のフォック演算子交換演算子(K)が、交換相関ポテンシャル(ν_XC)に変わっています。交換相関ポテンシャルに、ハートリー・フォック交換項に加えて電子相関の効果も含まれているため、DFTではより精度良く系の記述ができるようになっているそうです。

さらに、DFTはポスト・ハートリー・フォック法と比べて計算コストが低く抑えられています。ポストHF法では複数の電子配置を考えるので、精度を高めようと考慮する励起配置の数を増やしていくと、計算コストがどんどんと高くなっていきます。対してコーン・シャム法では「基底状態の電子配置のみを考え、かつ1電子方程式を解くのみで電子相関を取り込める」ので、計算が楽!ということのようです。

また、コーン・シャム近似では運動エネルギー汎関数の表現のためにコーン・シャム軌道を導入しました。 注意すべきなのは、この軌道は問題をとくために仮想的に導入された人工の系(コーン・シャムの補助系)における電子軌道であり、 「直接には物理的な意味をもっていない」という点です。従ってコーン・シャム方程式のエネルギー固有値も実際の系のものと見做すことはできないそうです。

ただし実際には、コーン・シャム方程式のエネルギー固有値は実験値と良く比較され、その根拠はJanakの定理と呼ばれるものにあるようです。より詳しくはこちらの記事をご参照ください。

「ハートリー・フォック方程式とコーン・シャム方程式似てるし、解き方も一緒だから、結果の意味も同じでしょ!」と思ってしまいましたが、そんな簡単な話ではなかったんですね!詳しくは理解できませんでしたがとりあえずコーン・シャムに基づいてエネルギーの議論に持ち込んでもOKそうなので安心(??)。

まとめ

以上、雰囲気量子化学入門(後編)でした。後編では、ハートリー・フォック行列の誤差(電子相関)から出発し、 その解決策として、拡張であるポスト・ハートリー・フォック法をみたあと、本命密度汎関数理論(DFT)の概要を辿りました。

イマイチ、出発点の動的電子相関と静的電子相関の具体的な意味を掴みきれていない気もしますが。。。

ぼんやりと聞いていたDFTが、波動関数ではなく電子密度に基づいたものであり、シュレーディンガー方程式ではなく等価な別の問題を解くものであったというのは全く理解していなかったので、びっくりしました。

別問題とはいいつつも、式の形や解くための変換の考え方はハートリー・フォックと共通する部分が多かったので、 前回の記事と合わせて全体の流れを辿ること何となく話を追うことができました。

中島隆人先生の「量子化学 ー分子軌道法の理解のためにー」、初学者にとって分かりやすい構成になっていてオススメですよー!

相変わらず、具体的な数式部分はさっぱり分かりませんでしたが、とりあえずB3LYP/6-31Gの意味がなんとなくわかったのでヨシ!!

今回もひどい誤解がたくさんあると思うのでご指摘いただければ幸いです。

参考資料

この記事は前編同様主に「量子化学 ー分子軌道法の理解のためにー 中島隆人 著」(裳華房)の内容を参考にしています。

上記書籍はハートリー・フォック法の問題点の指摘とその解決策の流れで密度汎関数理論が導入されており、流れを追いやすい一方で、 密度汎関数理論そのものについての記載量はあまり多くありません。 そこで、この記事の密度汎関数理論における記載はインターネット上に公開してくださっている以下の資料を参考にさせていただいています。 ただし、非専門家の私の独学での把握に基づいて記載しているので、誤りの責は全て私にあります。

密度汎関数理論の大まかな概要と用語についてはChem-Stationゼロさんの一連の記事がコンパクトで非常にわかりやすかったです。

より踏み込んだ内容、数式の展開、コーン・シャム法の意味するところなどについては@dc1394さんの一連の解説記事が非常に参考になりました。

また、汎関数についての具体的な例と説明、B3LYP/6-31Gについては以下の記事が参考になりました。

ハートリー・フォック法の誤差(電子相関)に関する記述は以下を参考にしています。動的電子相関・静的電子相関の具体的な意味合いについてはYouTube湯どうふさんの動画が大変参考になりました。

常田貴夫先生のこちらの本、Chem-Stationでも計算科学.comでもオススメされているのでとても気になるのですが、 懐が寂しいので買うべきか迷っています。。。

※ 書き終わってから気づきましたが、大阪大学 白井光雲先生の以下の資料は「原理の主張意味」に焦点をあてて解説してくださっていて、とてもわかりやすいかったです。もっと早く知りたかった。

*「第一原理と密度汎関数理論」 11 March 2005

おまけ(メモ)

前回、ハートーリー・フォック法では経験的なパラメーターを用いないので、HF法に基づく分子軌道法ab initio分子軌道法非経験的分子軌道法と呼ぶとかきました。 これに対して、煩雑な計算をすべて(1から)解くのではなく、一部、実験で得られたパラメーターをつかって計算する方法が半経験的分子軌道法です。 計算コストを下げる以外にも、電子相関の効果をパラメーターを通して取り込むといった点が期待できるそうです。 手法の例としてはAM1(Austin Model 1)法や、PM3(Parametic Methot 3)法といったものがあります。

  • v-表示可能性問題

DFTの基本定理「ホーヘンベルク・コーンの定理」の第1定理(存在定理)は「v-表示可能性の仮定」を用いています。これは「任意の電子密度 ρを与えると、必ず、それを基底状態の電子密度とするような外部ポテンシャル v が1つは存在する」というものです。 *12 この表現を使うと、存在定理は「基底状態波動関数と、v-表示可能である電子密度とは1対1の対応がある」と言い換えることもできるそうです。*13

問題は、「v-表示可能性であることの必要十分条件がわかっていない」ことで、このv-表示可能性問題の解決として「レヴィ(Levy)の制限付き探索」があります。

  • レヴィの制限付き探索

レヴィらは「v-表示可能な範囲」よりも広い「N-表示可能」な範囲に拡張して以下のように書き直しました。

f:id:magattaca:20210124135609p:plain

この密度汎関数Q[ρ] (Levy-Liebの密度汎関数)を使うとv-表示可能な領域の外でもホーヘンベルク・コーンの定理の第2定理(密度変分原理)が成り立ちます。これに従って変分探索を行い、エネルギーを最小化するような電子密度(ρ)を探せばv-表示可能性問題を回避できることになります。

もう少し具体的には、2段階の最小化工程となっています。

  1. まずρを固定して、そのρを与える波動関数Ψ_ρの中で「運動エネルギー(^T)と電子間相互作用(^V_ee)の期待値」を評価し、期待値を最小化するようなΨ_ρを探します。このΨ_ρをQ[ρ]と定義します。(1つ目の式)
  2. 次に、今度はρ固定せずに、2つ目の式のエネルギーを最小化するようなρを探します。

このプロセスが「レヴィの制限付き探索」と呼ばれるものです。

尚、N-表示可能は「N電子系において基底状態に対応する反対称化した波動関数(Ψ)を与える電子密度(ρ)」といったような意味で、以下のような簡単な条件をみたすどのようなρに対してもそのようなΨが存在することがわかっているそうです(電子密度のN表現可能性)。

f:id:magattaca:20210124134926p:plain

以上、よくわからなかったので本文に含めなかった事項の追加メモでした。

あと、いらすとやさんなんでもある。

*1: ハートリー・フォック極限のエネルギーで99.5%程度見積れるそうです。HF極限は基底関数を無限に増やしたHFR 計算によるエネルギーに対応するもののようです。

*2:以下の図で矢印の向きがごちゃごちゃしてしまってすみません。左半分について、「正確なエネルギーよりもハートリー・フォックエネルギーの方がエネルギーの値が高く、精度が上がるに連れてエネルギー値が低くなっていく」ように並べています。HF方程式を導くときに変分法を使ったので、「変分原理 :近似のエネルギー ε ≧ 真のエネルギー E 」を意識すれば、このような並び順となるはず。たぶん・・・(参考 wikipedia - Electronic correlation

*3:参考 Wikipedia - 電子相関

*4:電子相関の説明と図の作成に当たってはWebに公開されている理化学研究所 常田貴夫先生の講義資料とYouTube湯どうふ(化学の勉強)さんの説明動画を参考にさせていただいています

*5:波動関数は電子の存在確率に関係する(コペンハーゲン解釈)ということを踏まえている??

*6:価電子軌道(HOMO,LUMO)付近の分子軌道を用いて、その組の中で電子の励起で生じる電子配置を考慮(active space)します。CASSCFではこのactive spaceの全ての電子配置を考慮し、配置関数とする方法です。

*7:尚、多参照に対して単参照電子相関法はHF波動関数から出発して動的電子相関を考慮する相関法のことです。

*8:参考 できるだけ簡単に密度汎関数理論(Density Functional Theory, DFT)を説明してみる(前編)

*9:原子の全エネルギーを15%~50%程度低くみつもってしまうそう

*10: コーン・シャム軌道の項の表式(T_KS) 参照

*11:参考 「大学院講義 有機化学 Ⅰ.分子構造と反応・有機金属化学」 東京化学同人 第1版

*12:外部ポテンシャルは、物質中の電子に対しては原子核によるCoulombポテンシャルを意味するそうです。

*13:前記事のコメント欄でご紹介いただいた資料をご参照ください