230803
https://gyazo.com/d9ddd489379f4f1d2665f6502fa2436e
資料へのリンク
https://gyazo.com/099cba554c08c2e395a9fff063a8e07e
概要
日程： 2023 年 8 月 3 日（木）14時開始
会場： オンライン（Zoom）
スケジュール
table:schedule
開始 終了 時間 内容
14:00 14:05 5分 はじめに
14:05 14:30 25分 解説①：DFT計算の概要
14:30 14:55 25分 解説②：交換・相関汎関数をどう選ぶ？
14:55 15:20 25分 解説③：実践例と注意点
15:20 15:30 10分 質疑 / さいごに
15:30 終了
はじめに
学びの手引き
この資料について
資料については、講座が始まるまでに and/or 講座が終わってから、じっくりと読みながら学んでいただくことができるように準備しています
講座の終了後もできる限り公開を続けたいと思っているのですが、事情により、公開先を変更することがあるかも知れません
この講座の進め方
時間が限られてしまいますので、「ポイントのみを大まかに」という方針で説明することで、「DFT 計算に取り組むときの参考にしていただけそうな情報」をコンパクトにお伝えしたいと思っています
疑問に思われたこと、詳細が知りたいと思われたこと、その他のこと、どんなことでも、Zoom の「チャット」からコメントをお送りください
プレゼンテーションモード
資料を共有する時には、プレゼンテーションモードという機能を使っています
ページの右側にある を押して、Start presentationを選ぶと、プレゼンテーションモードになります
この講座の講師
山本典史 (Yamamoto, Norifumi)
https://res.cloudinary.com/dfhjad2jo/image/upload/v1672805236/yamnor_o9bsoi.png
Profile
専門は計算化学、コンピュータを使って分子を解析しています
化学の学びを身近にすることにも興味を持っています
連絡先
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください
Links
この講座の参考書・参考文献
参考書
参考文献
DFT計算の概要
復習：波動関数を用いる方法
原子や分子のふるまいを波動関数$ \Psi(\xi_1,\xi_2,...,\xi_n)を用いて解析する「おなじみの」方法を簡単に復習してみましょう。
Hamiltonian$ \hat{H}は、$ \Psiに対する演算子であり、$ n電子系では次のように表すことができます：
$ \hat{H} = - \frac{1}{2} \sum_i^n \nabla_i^2 - \sum_i^n v(\mathbf{r}_i) + \sum_{i>j}^n \frac{1}{r_{ij}}
ここで $ v(\mathbf{r_i})は、$ i番目の電子に外部（原子核）からはたらくポテンシャルなので、外部ポテンシャルと呼ぶことにします：
$ v(\mathbf{r}_i) = -\sum_a^{N_\mathbf{nuclei}} \frac{Z_a}{|\mathbf{R}_a - \mathbf{r}_i|}
https://gyazo.com/b208da7355872614f8cf35cc494ddd3e
Schrödinger 方程式 $ \hat{H} \Psi = E \Psiに基づくと、$ \Psiが決まれば、系のエネルギー$ Eを求めることができます：
$ E = \int \Psi^{*} \, \hat{H} \, \Psi \, \mathrm{d}\xi
たとえば、化学反応に伴う$ Eの変化が分かれば、反応速度などを理論的に予測することができるようになります。
DFTとは？
密度汎関数法（DFT）は、波動関数$ \Psi(\xi_1,\xi_2,...,\xi_n)の代わりに電子密度$ \rho(\mathbf{r})を用いて、原子や分子のふるまいを解析する方法です。
電子を、粒子でも波動でもなく、密度雲（密度分布）という概念で捉える。
$ \Psi(\xi_1,\xi_2,...,\xi_n)は、$ n電子系の場合、$ 3n個の変数$ \xi=x_i,y_i,z_i \;(i=1,2,...,n)を含んでいます。
➡ 電子数$ nが多くなると、計算が大変！
一方で$ \rho(\mathbf{r})の場合、変数は$ \mathbf{r}(x,y,z)の 3 個のみになります。
➡ $ \Psiの代わりに$ \rhoを使えば、変数が$ 3nから 3 に減少して、$ nが多くなっても大丈夫ですね！
まとめると...
https://gyazo.com/f3b401fbb48d795614ae3ef5ac4b808c
ただし、実際にDFTを活用する場合、 1 電子波動関数（軌道）を使った方法（Kohn-Sham 法）に頼らざるを得ないので、DFT 計算に必要な変数は結局、 $ 3n個に戻ってしまいます...
DFTの基本原理①
Hohenberg-Kohnの第一定理によると、$ \rhoが与えられれば、基底状態における原子や分子の全電子エネルギー$ Eを求めることができます。
https://gyazo.com/4a5f58367034e286252aadfd91cafc67
この場合、$ Eは$ \rhoの汎関数であるといい、次のように表します：
$ E = E[\rho]
汎関数って何？
関数は、数値を入力すると、数値を出力するもの。
たとえば$ \rho(\mathbf{r})は、$ \mathbf{r}を入力として、その座標における電子密度の値を出力する関数です。
汎関数は、関数を入力して、数値を出力するもの。
$ E[\rho(\mathbf{r})] は、関数である$ \rho(\mathbf{r}) を入力として、その電子密度に対応するエネルギーの値を出力する汎関数です。
DFTの基本原理②
Hohenberg-Kohnの第二定理によると、ある系の基底状態における真のエネルギー$ E_0は、ある試行的な電子密度$ \rho_\mathrm{trial}に対するエネルギーに対して、次のような関係が成り立ちます：
$ E_0 \leqq E[\rho_\mathrm{trial}]
➡ なるべく$ Eが小さくなるような$ \rho_\mathrm{trial}を探せばよい！
DFTの課題
DFTの基本原理にしたがうと、基底状態における原子や分子のエネルギー$ Eは、適切な$ \rho(\mathbf{r})から求めることができるはず、ですが...
$ E = E_v[\rho(\mathbf{r})] = T[\rho] + E_\mathrm{ne}[\rho] + J[\rho] + E_\mathrm{x}[\rho]
$ T[\rho] ：運動エネルギー項
➡ 具体的なかたちが分かりません
$ E_\mathrm{ne}[\rho] ：核-電子 Coulomb 相互作用項
➡$ E_\mathrm{ne}[\rho] = \int v(\mathbf{r})\rho(\mathbf{r}) \;\mathrm{d}\mathbf{r}
$ J[\rho] ：古典的な電子-電子「 Coulomb 反発」項
➡$ J[\rho] = \frac{1}{2} \iint \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}
$ E_\mathrm{x}[\rho] ：非古典的な電子-電子「交換」項
➡ 具体的なかたちが分かりません
Thomas-Fermi-Dirac（TFD）モデル
DFT の課題は...
$ T[\rho] と$ E_\mathrm{x}[\rho] の具体的なかたちが分からない
$ \rho({\mathbf{r})}のふるまいに「なんらかの制限条件」を付ければ、近似式が得られるかも。
Thomas と Fermi は、ある系が局所的には一様密度の電子ガス（均一電子ガス）としてふるまうと近似した場合に対して、運動エネルギー汎関数を次のように導出しました：
$ T^\mathrm{TF}[\rho] = \frac{3}{10}(3\pi^2)^{2/3} \int \rho^{5/3}(\mathbf{r}) \; \mathrm{d}\mathbf{r}
Dirac は、均一電子ガスにおける交換汎関数として、次の形式を提案しました：
$ E_\mathrm{x}^\mathrm{LDA}[\rho] = \frac{3}{4}\left( \frac{3}{\pi} \right)^{1/3} \int \rho^{4/3}(\mathbf{r}) \; \mathrm{d}\mathbf{r}
$ T[\rho] と$ E_\mathrm{x}[\rho] の具体的なかたちが分かって一件落着？？
均一電子ガスモデルは、孤立原子や分子に対して、適切な記述だとは言えません...
化学結合を全く考慮できません...
しかし、$ v(\mathbf{r})が空間的に一定ではないリアルな系にDFTを展開するときのベースになっています。
とくに、$ T[\rho] の記述が悪い...
運動エネルギーを、一般的に、10% 程度も過小評価する（らしい）。
Kohn-Sham 法
TFD モデルの大きな欠点は...
運動エネルギー汎関数$ T[\rho] の記述が悪い
この問題に対して Kohn と Sham は...
電子密度$ \rho 使わず、1 電子波動関数（軌道）$ \phi を使って$ T[\rho] を表現するアイデアを提唱しました。
このアイデアを大まかに説明すると...
有効外部ポテンシャル $ v_s(\mathbf{r})のもとにある $ n個の相互作用がない仮想の電子系（Kohn-Sham 系）を考えます。
ただし$ v_s(\mathbf{r})は、相互作用を考慮したリアルな系で求められるものと等しい（仮想的な分子系）とします：
https://gyazo.com/9d050d7261c04608e767abdd31023288
Kohn-Sham 方程式
Kohn-Sham 系に対する Schrödinger 方程式（Kohn-Sham 方程式）は、次のようになります：
$ \hat{h}_i^\mathrm{KS} \phi_i^\mathrm{KS} = \epsilon_i^\mathrm{KS} \phi_i^\mathrm{KS}
$ \hat{h}_i^\mathrm{KS} = -\frac{1}{2} \nabla_i^2 + v_s(\mathbf{r}_i)
Kohn と Sham は、KS 軌道 $ \phi_i^\mathrm{KS}を用いて、運動エネルギー汎関数を次のように表しました：
$ T_s[\rho] = -\frac{1}{2}\sum_i \langle \phi_i^\mathrm{KS} \, | \, \nabla_i^2 \, | \, \phi_i^\mathrm{KS} \rangle
エネルギー汎関数は、$ T_s[\rho] を使うと、
$ E_v[\rho] = T_s[\rho] + \frac{1}{2} \iint \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}+ E_\mathrm{xc}[\rho] + \int v(\mathbf{r})\rho(\mathbf{r}) \;\mathrm{d}\mathbf{r}
最後に残ったのは、交換・相関汎関数$ E_\mathrm{xc}[\rho]
$ E_\mathrm{xc}[\rho] を厳密に求める方法は分かっていないので、さまざまな種類の近似汎関数が提案されています。
交換・相関汎関数をどう選ぶ？
交換・相関汎関数とは？
DFT の計算精度と信頼性は、近似汎関数$ E_\mathrm{xc}[\rho] の種類に大きく依存します。
$ E_\mathrm{xc}[\rho] を、交換項 $ E_\mathrm{x}[\rho] と相関項$ E_\mathrm{c}[\rho] に分けて表すと...
$ E_\mathrm{xc}[\rho] = E_\mathrm{x}[\rho] + E_\mathrm{c}[\rho]
DFT の黎明期には、まずは$ E_\mathrm{x}[\rho] と$ E_\mathrm{c}[\rho] を個別に開発し、あとからそれらを適切に組み合わせていました。
現在では、最初から２つを組み合わせた$ E_\mathrm{xc}[\rho] の開発に取り組むことが主流のようです。
交換・相関汎関数の分類
https://gyazo.com/f9a0046baad0a2b25cf7d682c8f9adbb
局所スピン密度近似（LSDA）型
電子密度 $ \rho のみで表現された汎関数です。
一般化勾配近似（GGA）型
LSDA を電子密度勾配 $ \nabla \rho で補正した汎関数です。
meta GGA 型
GGA を運動エネルギー密度 $ \tauで補正した汎関数です。
hybrid 型
Hartree-Fock 交換積分$ E_\mathrm{x}^\mathrm{HF}を一定量混合した汎関数です。
LSDA 交換汎関数
密度が局所的に均一電子ガス（UEG）として扱えると仮定します。
https://gyazo.com/afda28afe877bafc7c846f9c98c45c45
TF モデルで説明したように、UEG に対する交換汎関数として、Dirac が次のものを提唱していました：
$ E_\mathrm{x}^\mathrm{LDA}[\rho] = \int e_{\mathrm{x}}^\mathrm{UEG}(\rho) \; \mathrm{d}\mathbf{r}
$ e_\mathrm{x}^\mathrm{UEG}(\rho) = -\frac{3}{4}\left( \frac{3}{\pi} \right)^{1/3} \rho^{4/3}
この汎関数は、空間の各点で（つまり局所的に）、その場所の電子密度$ \rho(\mathbf{r})だけの関数になっていることから、局所密度近似（LDA）と呼ばれます。
https://gyazo.com/fbc31958ac35ce0c52c1334a2b769b49
上記について、スピン偏極した系に対応させたものが LSDA 交換汎関数です。
$ E_\mathrm{x}^\mathrm{LSDA}[\rho] = \sum_\sigma^{\alpha,\beta} \int e_{\mathrm{x},\sigma}^\mathrm{UEG}(\rho) \; \mathrm{d}\mathbf{r}
$ e_{\mathrm{x},\sigma}^\mathrm{UEG}(\rho) = -\frac{3}{2}\left( \frac{3}{4\pi} \right)^{1/3} \rho^{4/3}
LSDA は、UEG というモデル系に対する厳密な DFT になります。
金属のように、電子密度が緩やかに変化するような系を扱う場合には、広く使われているようです。
分子系に対しては、交換エネルギーを 10% 程度も過小評価してしまいます。
GGA 交換汎関数
GGA は、原子や分子のような電子密度が均一ではない系を扱うために、LSDAに対する補正として、$ \rho の一階微分$ \nabla \rho を変数として取り込んだ汎関数です。
$ E_\mathrm{x}^\mathrm{GGA}[\rho,\nabla\rho] = \sum_\sigma^{\alpha,\beta} \int e_{\mathrm{x},\sigma}^\mathrm{UEG}(\rho) \; g_{\mathrm{x},\sigma}^\mathrm{GGA}(\rho,\nabla\rho) \; \mathrm{d}\mathbf{r}
$ g_{\mathrm{x},\sigma}^\mathrm{GGA}は、不均一性補正係数（ICF）と呼ばれます。
さまざまな$ g_{\mathrm{x},\sigma}^\mathrm{GGA}の形式が提唱されています。
例：OPTX, PBE, B86, B88, PW86, PW91, ...
無次元の勾配変数 $ s=\frac{1}{2(2\pi^2)^{1/3}}\frac{|\nabla\rho|}{\rho^{4/3}} に対する$ g_{\mathrm{x},\sigma}^\mathrm{GGA}の変化を図示してみると：
補足：LSDAは、常に$ g_x=1となります。
https://gyazo.com/b310b915c38b3f9056df2ed76c4af420
meta GGA 交換汎関数
meta GGA は、GGA を拡張して、電子密度の高次導関数（$ \nabla^2 \rho）に依存するようにしたものです。
運動エネルギー密度$ \tauは、電子密度のラプラシアン（$ \nabla^2 \rho）と本質的に同じ情報を持ちます。
meta GGA では、一般に、数値的に安定な$ \tauの方を変数に用います。
$ E_\mathrm{x}^\mathrm{mGGA}[\rho,\nabla\rho,\tau] = \sum_\sigma^{\alpha,\beta} \int e_{\mathrm{x},\sigma}^\mathrm{UEG}(\rho) \; g_{\mathrm{x},\sigma}^\mathrm{mGGA}(\rho,\nabla\rho,\tau) \; \mathrm{d}\mathbf{r}
さまざまな$ g_{\mathrm{x},\sigma}^\mathrm{mGGA}の形式が提唱されています。
例：PTZB, TPSS, MS2, MVS, SCAN, ...
無次元の勾配変数 $ s=\frac{1}{2(2\pi^2)^{1/3}}\frac{|\nabla\rho|}{\rho^{4/3}} に対する$ g_{\mathrm{x},\sigma}^\mathrm{mGGA}の変化を図示してみると：
補足：LSDAは、常に$ g_x=1となります。
$ \alpha=0：共有結合に対して
$ \alpha=\infin：非共有結合（分子間結合, etc）に対して
https://gyazo.com/c9f8aca6f9619c848135a9bf0a40580b
hybrid 交換・相関汎関数
hybrid 型は、LSDA / GGA / mGGA などの交換汎関数に対して、一定の割合で Hartree-Fock 交換積分$ E_\mathrm{x}^\mathrm{HF}[\phi] を混合した汎関数です。
$ E_\mathrm{x}^\mathrm{hybrid}[\rho,\phi] = \int_0^1 E_\mathrm{x}^\lambda \; \mathrm{d}\lambda = E_\mathrm{x}^\mathrm{DFT}[\rho] + \lambda (E_\mathrm{x}^\mathrm{HF}[\phi] - E_\mathrm{x}^\mathrm{GGA}[\rho])
上記は、断熱接続式（ACF）と呼ばれています。
HF 交換積分 $ E_\mathrm{x}^\mathrm{HF}は、KS 軌道$ \phi_i^\mathrm{KS}を用いて、次のように求めることができます：
$ E_\mathrm{x}^\mathrm{HF}[\phi^\mathrm{KS}] = -\frac{1}{4} \sum_{i,j}^n \langle \phi_i^\mathrm{KS}(\mathbf{1}) \phi_j^\mathrm{KS}(\mathbf{2}) \; | \; 1/r_{ij}\; | \; \phi_j^\mathrm{KS}(\mathbf{1}) \phi_i^\mathrm{KS}(\mathbf{2}) \rangle
$ \lambdaは、電子-電子相互作用の程度を決めるパラメータです。
$ \lambda = 0の極限では、独立電子モデル（KS モデル）に基づく交換汎関数となります。
たとえば、代表的な交換・相関汎関数である B3LYP のかたちは、次のように表すことができます：
$ E_\mathrm{xc}^\mathrm{B3LYP} = (1-a) E_\mathrm{x}^\mathrm{LSDA} + a E_\mathrm{x}^\mathrm{HF} + b E_\mathrm{x}^\mathrm{B88} + (1-c) E_\mathrm{c}^\mathrm{LSDA} + c E_\mathrm{c}^\mathrm{LYP}
$ a = 0.20, \; b=0.72, \; c=0.81
➡ ３個のパラメータは、数十種類の原子＆小さな分子について、その物性値を再現するようにフィッティング
非 hybrid 型（local 型）と hybrid 型の汎関数ペア（例：BLYP vs B3LYP）について、計算される物性値の平均二乗偏差（RMSD）を比較してみると...
NCD：分子間相互作用エネルギー
ID：異性化エネルギー
BH：反応障壁の高さ
https://gyazo.com/eecd63c52ff0735df7ece58b720e7b39
$ E_\mathrm{x}^\mathrm{HF}を適度に混合させることで計算精度が向上しています。
改善の一部は、$ E_\mathrm{x}^\mathrm{HF}によって、自己相互作用誤差（次で説明）が減少するためだと考えられます。
自己相互作用誤差とは？
水素原子（系に含まれる電子が１個のみ）の場合...
Coulomb エネルギー $ J[\rho] は、系に含まれる１個の電子自身の自己相互作用になります。
自己相互作用は、対応する交換エネルギー$ E_\mathrm{xc} [\rho] によって打ち消されるはずです。
つまり、$ J[\rho] + E_\mathrm{xc} [\rho] = 0 である必要があります。
HF 法の場合、完全基底系（CBS）で計算すると...
Coulomb エネルギー：$ J_{ii}[\phi] = 0.3125 Hartree
交換エネルギー：$ K_{ii}[\phi] = 0.3125 Hartree
HF 法は、常に $ \sum_i^n (J_{ii}[\phi] - K_{ii}[\phi])=0 であり、 自己相互作用誤差を含んでいません。
DFT の場合、一般に...
$ \Delta E^\mathrm{SIE} = \sum_i ( J[\rho(\mathbf{r}_i)] + E_\mathrm{xc} [\rho(\mathbf{r}_i)]) \neq 0 であり、自己相互作用誤差を含みます。
➡ 電子が不自然に局在する傾向があり、エネルギーの過大評価が起こってしまいます。
hybrid 型の交換・相関汎関数を使うと...
HF 交換積分$ E_\mathrm{x}^\mathrm{HF}を含むので、$ \Delta E^\mathrm{SIE} が減少し、計算精度が向上します。
交換・相関汎関数の階層的な分類
Perdew によって提案された汎関数の階層的な分類です。
https://gyazo.com/70fdb4eec4ffaa4eee9e5d8643beecaf
はしごを１段上がるごとに計算精度が向上すると期待されますが、そうなる保証はありません。
一般的によく利用されているのは、Level 2 にある GGA 汎関数、Level 4 にある hybrid 型です。
交換・相関汎関数の妥当性
交換・相関汎関数には、それぞれ、長所と短所があります。
現状、対象とする系や解析しようとする内容に応じて、「最適」な汎関数を適切に選ぶ必要があります。
https://gyazo.com/8c6b2ffb60f76b70b816a1c6c3dca8e7
汎関数に含まれる半経験的なパラメータ
B97 系や Mx 系など、最近開発されている汎関数の多くは、実験データとよく一致するようにパラメータをフィッティングして、計算精度を向上させようとしています。
これから計算しようとする系が、パラメータを決めるときに使った系と似たものである場合、高い計算精度を与える傾向にあります。
一方で、パラメータを決めるときに利用するデータは第３周期までの原子（H 〜 Ar）から構成される分子系であることが多く、遷移金属を含む場合には、奇妙なふるまいを示すことがあります。
https://gyazo.com/35c79c38d998857dc5a969e60389113a
いくつかの汎関数には、40〜60 個という、膨大な数のフィッティングパラメータが含まれています。
半経験的なパラメータを使った汎関数は、「洗練された汎関数が開発されるまでの過渡的な汎関数である」（Tsuneda (2012)）という意見もあります。 交換・相関汎関数に対する補正
長距離補正
交換汎関数で考慮されていない長距離電子間の交換相互作用に対する補正です。
分散力補正
相関汎関数で考慮されていない分散力に対する補正です。
長距離補正
交換汎関数は、一般に、長距離電子間の交換相互作用が考慮されていません。
理由①：電子密度 $ \rho(\mathbf{r}) は、２電子間のあらわな相互作用に関する情報は含んでいないから。
理由②：一般的な交換汎関数は、１電子座標のみに関する積分で計算されるから。
例：$ E_\mathrm{x}^\mathrm{GGA}[\rho] = \sum_\sigma^{\alpha,\beta} \int e_{\mathrm{x},\sigma}^\mathrm{UEG}(\rho) \; g_{\mathrm{x},\sigma}^\mathrm{GGA}(\rho,\nabla\rho) \; \mathrm{d}\mathbf{r}
一方で、HF 交換積分 $ E_\mathrm{x}^\mathrm{HF}は、２電子積分のなかで２電子間の相互作用をあらわに計算するため、長距離相互作用が自然なかたちで考慮されます。
$ E_\mathrm{x}^\mathrm{HF}[\phi] = -\frac{1}{2} \sum_\sigma^{\alpha,\beta} \sum_{i,j}^n \langle \phi_{i,\sigma}^\mathrm{KS}(\mathbf{1}) \phi_{j,\sigma}^\mathrm{KS}(\mathbf{2}) \; | \; \frac{1}{r_{ij}} \; | \; \phi_{j,\sigma}^\mathrm{KS}(\mathbf{1}) \phi_{i,\sigma}^\mathrm{KS}(\mathbf{2}) \rangle
長距離相互作用誤差が問題になるのは...
van der Waals 相互作用
電子励起スペクトルや光学応答性（特に、電荷移動状態や Rydberg 状態を含むとき）
そこで...
短距離（sr）部分を $ E_\mathrm{x}^\mathrm{DFT}[\rho] 、長距離（lr）部分を $ E_\mathrm{x}^\mathrm{HF}[\phi] で表すように領域分割する長距離補正（LC）が提案されました。
hybrid 型に対して長距離補正をおこなった汎関数は、領域分割 hybrid（RSH）型といい、次のように表します：
$ E_\mathrm{xc}^\mathrm{RSH} = c_\mathrm{x,sr} E_\mathrm{x,sl}^\mathrm{HF} + c_\mathrm{x,lr} E_\mathrm{x,lr}^\mathrm{HF} + (1-c_\mathrm{x,sr}) E_\mathrm{x,sl}^\mathrm{DFT} + (1-c_\mathrm{x,lr}) E_\mathrm{x,lr}^\mathrm{DFT}
領域分割のふるまいを図示すると...
補足：y 軸が 1 のときに、HF 交換積分と等しくなります。
https://gyazo.com/cb1edc399ac69fd3458779b3591c4a73
上記の図の補足
領域分割は、誤差関数 $ \mathrm{erf} を用いて、二電子演算子 $ 1/r_{12}を次のように表すことでおこなっています：
$ \frac{1}{r_{12}}=\frac{1-\mathrm{erf}(\omega r_{12})}{r_{12}} + \frac{\mathrm{erf}(\omega r_{12})}{r_{12}}
分散力補正
標準的な DFT の弱点のひとつは、分散力（van der Waals 相互作用）の記述ができないことです。
➡ このために、分子集合体における分子間相互作用を過小評価する傾向にあります。
Grimme は、分子力場などで用いられるように経験的なパラメータ$ C_nを用いて、各原子対に対する引力的なエネルギー項 $ R^{-n} を付加する分散力補正を提案しました：
$ \Delta E_\mathrm{disp}=-\sum_{n=6,(8,10)}s_n\sum_{AB}^\mathrm{atoms}\frac{C_n^{AB}}{R_n^{AB}}f_\mathrm{damp}(R_{AB})
この方法を用いることによる追加の計算コストはほぼゼロです。
その他に、van Voorhis-Vydron（VV10）核という分散力カーネル $ \Phi(\mathbf{r},\mathbf{r}')を用いて、実際の電子密度に直接依存させる分散力補正も提案されています：
$ \Delta E_\mathrm{disp}=\frac{1}{2}\int\rho(\mathbf{r}) \, \Phi(\mathbf{r},\mathbf{r}') \, \rho(\mathbf{r'}) \, \mathrm{d}\mathbf{r} \, \mathrm{d}\mathbf{r}'
Grimme の分散力補正のあり・なし（例：BLYP と BLYP-D2, BLYP-D3）について、計算される物性値の平均二乗偏差（RMSD）を比較してみると...
NCED：分子間相互作用エネルギー
IE：異性化エネルギー
https://gyazo.com/260c089bbd967b6b9b70a2d5f55d930e
Grimme の分散力補正を加えることで、コストを追加することなく、計算精度を向上させることができます。
交換・相関汎関数の例
代表的な200個の交換・相関汎関数をまとめると...
Local：非 hybrid 型
Disp：分散力補正あり
下線：長距離補正あり
https://gyazo.com/7f3a9a074aae8c94477556d2bb3c3da5
交換・相関汎関数のベンチマーク
200 個の汎関数について、ベンチマーク計算をしてみると...
汎関数名の右にある数字は、計算精度のランキング
左にあるラベルは、評価した物性値：
NCED / NCED / NCD：分子間相互作用エネルギー
IE / ID：異性化エネルギー
TCE / TCD：熱化学量
BH：反応障壁の高さ
https://gyazo.com/42c648b2b9ebf15b3674168b8b6de3fc
GGA から meta GGA に変えても、多くの場合、計算精度はそれほど良くはならないようです。
非 hybrid 型（local 型）から hybrid 型に変えると、ほとんどの場合、計算精度が向上します。
計算例①：水分子の構造最適化
https://gyazo.com/c2e381a0aff6f0bfb93bac93b78c270f
計算例②：水分子での振動数解析
https://gyazo.com/288677666289ced8a6fec15061e4bcb8
対称・逆対称伸縮振動の誤差が大きいのは、分子振動の非調和性が原因です。
➡ 電子状態計算の計算レベルを上げるだけではなく、分子振動の非調和性を考慮する必要があります。
実践例と注意点
DFT を用いて励起状態を解析できる？
DFT は、基底状態の計算に関しては、HF 法と同等のコストで、post HF 精度の結果が得られる可能性があります。
DFT は、原子や分子の励起状態に関しても、時間依存（TD）の問題として扱うように拡張（TD-DFT）できます。
TD-DFTは、化学発光 / 光触媒 / 太陽電池 / 光合成など、さまざまな分野で活用されています。
ただし、TD-DFT を用いて励起状態を解析するときには、注意が必要な場合もあります。
今回は、TD-DFT の概要を説明したあと、具体的な実践例と注意点について簡単に紹介します。
TD-DFTとは？
時間依存 Kohn-Sham 方程式は、時間依存 Schrödinger 方程式と同様に、次のように表すことができます：
$ \hat{h}_i^\mathrm{KS} \phi_i^\mathrm{KS} = i\frac{\partial}{\partial t} \phi_i^\mathrm{KS}
ここで...
$ \hat{h}_i^\mathrm{KS} = -\frac{1}{2} \nabla_i^2 + v_\mathrm{ext}(\mathbf{r}_i,t) + J(\mathbf{r}_i,t) + v_\mathrm{xc}(\mathbf{r}_i,t)
Coulomb ポテンシャル：$ J(\mathbf{r}_i,t) = \int\frac{\rho(\mathbf{r}',t)}{|\mathbf{r} - \mathbf{r}'|} \; \mathrm{d}\mathbf{r}'
交換相関ポテンシャル：$ v_\mathrm{xc}(\mathbf{r}_i,t) = \frac{\delta E_\mathrm{ex}[\rho(\mathbf{r}_i,t)]}{\delta\rho(\mathbf{r}_i,t)}
線形応答理論に基づいて、時間依存 KS 方程式を解くと...
$ \def\mat#1{\begin{bmatrix}#1\end{bmatrix}} \mat{\mathbf{A}&\mathbf{B}\\\mathbf{B}^{*}&\mathbf{A}^{*}}\mat{\mathbf{X}\\\mathbf{Y}}=\omega\mat{\mathbf{1}&\mathbf{0}\\\mathbf{0}&-\mathbf{1}}\mat{\mathbf{X}\\\mathbf{Y}}
$ \pm \omegaの値が、励起 ＆ 脱励起のエネルギーに対応します。
行列要素は...
$ A_{ij}^{ab} = \delta_{ij}\delta_{ab}(\epsilon_a - \epsilon_i) + \langle ij | ab\rangle + \langle ij \, | \, f_\mathrm{xc} \, | \, ab\rangle
$ B_{ij}^{ab} = \langle ij | ab\rangle + \langle ij \, | \, f_\mathrm{xc} \, | \, ab\rangle
行列 $ \mathbf{B}を無視して上記の方程式を解く扱いのことを、Tamm-Dancoff 近似（TDA）といいます：
$ \mathbf{AX} = \omega \mathbf{X}
➡ 遷移モーメントの定量性が悪くなる可能性があります。
➡ 一方で、励起エネルギーの計算精度が改善されることがあります。
$ f_\mathrm{xc}は、交換・相関カーネルといい、次のように表します：
$ f_\mathrm{xc}(\mathbf{r},t,\mathbf{r}',t') = \frac{\delta v_\mathrm{xc}(\mathbf{r},t)}{\delta \rho(\mathbf{r}',t)}
TD-DFT のベンチマーク（TDA など）
GINV：meta GGA 汎関数 $ E_\mathrm{xc}^{mGGA}[\rho,\nabla\rho,\tau] の運動エネルギー密度$ \tau(\mathbf{r},t)のゲージ依存性を補正した結果
https://gyazo.com/419563c10e4a060c1fffa55a97c10b35
三重項状態の計算では、TDA によって計算精度が改善しています。
TDA はゲージ不変性を破ってしまうので、GINV と同時に採用することはできません。
GINV による改善はわずかなので、通常は、TDA の方をオススメします。
TD-DFT のベンチマーク（汎関数依存性）
https://gyazo.com/53e8a78345710b4d969c72d291f8d929
領域分割 hybrid 型（長距離補正）を用いることで、計算精度が改善しています。
蛍光（fluorescence）の結果に関しては、ベンチマークに用いたデータが少ないので、信頼性に欠けるかもとのこと。
TD-DFT の課題
一般に、TD-DFT で交換・相関カーネル$ f_\mathrm{ex}(\mathbf{r},t,\mathbf{r}',t')を扱う場合、その時間依存を無視します（断熱近似）：
$ f_\mathrm{xc}^\mathrm{adia}(\mathbf{r},\mathbf{r}') = \frac{\delta v_\mathrm{xc}(\mathbf{r})}{\delta \rho(\mathbf{r}')} = \frac{\delta^2 E_\mathrm{xc}[\rho(\mathbf{r})]}{\delta \rho(\mathbf{r})\delta \rho(\mathbf{r}')}
$ f_\mathrm{ex}^\mathrm{adia}を用いる場合、１電子励起のみを考慮した描像になることに注意が必要です。
さまざまな光化学現象を記述するときに必要な２電子励起の寄与を考慮できません。
断熱近似は、多くの場合、妥当な結果を与えるのですが、いくつかのケースでは問題もあります。
共役系の低い励起状態の計算精度が悪くなることがあります。
CASPT2：6.27 eV
TD-DFT（B3LYP）：7.02 eV
基底状態と励起状態が交差する円錐交差を適切に扱うことができません。
円錐交差とは？
分子の立体構造が変化しながら励起状態と基底状態の間のエネルギー差が段々と小さくなり、最終的に、そのエネルギー差がゼロになってしまう場合があります。
この「電子状態間のエネルギー差がゼロとなって互いに交差する」ことを円錐交差といいます。
https://gyazo.com/3c2e37c2244f77ffce0c458640152629
円錐交差する地点に辿り着いたとき、分子は、光を放出しないで基底状態に戻ります。
励起状態となった分子が「光の放出をともなわずに基底状態に戻る」ことを無輻射遷移 or 非断熱遷移といいます。
非断熱遷移は、さまざまな機能性分子が多彩で巧みな機能を発現するための源です！
分子の光化学反応を解析するときに、円錐交差の地点を探索し、その立体構造を調べたり、非断熱遷移の起こりやすさを調べたいと思うことがしばしばあります。
円錐交差点の特徴
円錐交差点は、$ N原子系の場合、自由度 $ M=3N-6に対して、$ M-2次元空間にあるポテンシャルエネルギー曲面間の縮退点です。
円錐交差点は、一次元の「点」ではなくて、多次元の「縫い目（シーム）」
円錐交差点において、系のエネルギー縮退を解く２つのベクトル $ \mathbf{x}_1, \mathbf{x}_2を、分岐ベクトルといいます。
https://gyazo.com/42e44616bb7ef60e7d2e5fc3b076554a
エチレンの円錐交差点における分岐ベクトル
たとえば、エチレン（$ \mathsf{CH_2=CH_2}）は、二重結合の回転にともなって$ \mathsf{S_0}状態と$ \mathsf{S_1}状態が近接し、円錐交差点に至ります。
このとき...
分岐ベクトル$ \mathbf{x}_1：二重結合の回転方向
分岐ベクトル$ \mathbf{x}_2：ジラジカル構造 ↔ イオン構造の変化
https://gyazo.com/a6b2e9a5d8a60fc3c011593383061bee
TD-DFT はエチレンの円錐交差点を正しく記述できない
エチレンの円錐交差点における分岐ベクトル $ \mathbf{x}_1, \mathbf{x}_2に沿ったエネルギー変化を見ると...
SSR：State-interection State-averaged spin-Restricted ensemble KS
SF：Spin-Flip TD-DFT
TD：TD-DFT
https://gyazo.com/2c5c1e8a8dd3fed2e56ec86e7f6c4a19
TD-DFT（断熱近似）の場合、円錐交差の構造を適切に記述できていません。
円錐交差点を適切に記述するには
DFT
スピン反転（spin-flip）TD-DFT
アンサンブル（ensemble）DFT
WFT
CASSCF などの多配置 SCF 法
SF-TD-DFT を用いた円錐交差の解析事例
まとめると...
円錐交差点は、ちょっとマニアックな話題に聞こえますが、光機能性材料の研究開発で意外に重要です。
凝集誘起発光という光化学現象のふるまいにも、円錐交差点が重要な役割を果たしています。
spin-flip TD-DFT を使うと、円錐交差点近傍の電子状態を記述することができます。
ただし spin-flip 法には、「スピン混入」という悩ましい問題があって、注意が必要です。
量子コンピュータを用いた円錐交差の解析事例
SA-OO-VQE アルゴリズムを用いた解析例：
https://gyazo.com/4a606135ec1b3d27573d87374f845728
質疑
さいごに
共同研究・技術相談
計算化学（分子動力学・量子化学・機械学習）に関して、共同研究や技術相談をお受けしています。
興味・関心を持っていただけましたら、山本までお問い合わせください。
お問い合わせ
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください