230824
https://gyazo.com/cfe1a52590eb89e1c5a9d41026c0f805
資料へのリンク
https://gyazo.com/9742a9f8fbde4e2db443a669d6e4aa36
概要
日程： 2023 年 8 月 24 日（木）
会場： オンライン（Zoom）
目次
はじめに
量子化学計算の概要
量子化学計算の代表的なソフトウェア
量子化学計算の代表的な計算方法
量子化学計算の代表的な基底関数
量子化学計算に基づく構造最適化
量子化学計算に基づく振動解析
量子化学計算に基づく化学反応の解析
量子化学計算に基づく電子励起状態の解析
質疑応答
さいごに
はじめに
学びの手引き
この資料について
資料は、講座が始まるまでに and/or 講座が終わってから、じっくりと読みながら学んでいただくことができるように準備しています
講座の終了後もできる限り公開を続けたいと思っていますが、事情により、公開先を変更することがあるかも知れません
この講座の進め方
時間が限られるので、「ポイントのみを大まかに」という方針で説明し、「量子化学計算をはじめようとするときに参考にしていただけそうな情報」をコンパクトにお伝えしたいと思っています
星印（★）の項目は応用的な内容です
「はじめて量子化学計算に取り組もうとするとき」には、スキップしておいても良いように思います
疑問に思われたこと、詳細が知りたいと思われたこと、その他のこと、どんなことでも、Zoom の「Q&A ツール」からコメントをお送りください
代表的なツールやアプリの「実演」をおこないながら、「量子化学計算を実際におこなうときのイメージを掴んでいただくこと」も大切にしたいと思っています
プレゼンテーションモード
資料を共有する時には、プレゼンテーションモードという機能を使っています
ページの右側にある を押して、Start presentationを選ぶと、プレゼンテーションモードになります
この講座の講師
山本典史 (Yamamoto, Norifumi)
https://res.cloudinary.com/dfhjad2jo/image/upload/v1672805236/yamnor_o9bsoi.png
Profile
専門はコンピュータ化学、コンピュータを使って分子を解析しています
化学の学びを身近にすることにも興味を持っています
連絡先
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください
Links
この講座の到達目標
量子化学計算に取り組むときに必要となる基礎的な知識と基本的な技術について、説明できる
量子化学の知識や技術に基づいて、おもに有機化合物の化学的性質を解析するための指針について、説明できる
代表的な量子化学計算ソフトウェアである Gaussian と GAMESS の特徴と基本的な計算手順について、説明できる
代表的な量子化学計算ソフトウェアのインターフェースとして活用できる Web アプリ WebMO の特徴と基本的な使い方について、説明できる
この講座のタイムスケジュール
table:schedule
開始 終了 時間 内容
10:30 12:00 90分 講義①
12:00 12:45 45分 昼食
12:45 14:15 90分 講義②
14:15 14:30 15分 休憩
14:30 16:00 90分 講義③
16:00 16:30 30分 質疑
16:30 終了
講義①
はじめに
量子化学計算の概要
量子化学計算の代表的なソフトウェア
講義②
量子化学計算の代表的な計算方法
量子化学計算の代表的な基底関数
講義③
量子化学計算に基づく構造最適化
量子化学計算に基づく振動解析
量子化学計算に基づく化学反応の解析
量子化学計算に基づく電子励起状態の解析
質疑
質疑応答
さいごに
この講座の参考書
はじめて量子化学計算に取り組む人向けとして、読みやすい入門書です
Gaussian に準拠した計算化学の入門書です
一般の書店では販売されていません
Gaussian の国内代理店などで購入できます
量子化学の基礎を学ぶ
量子化学計算の基礎となる「量子論」から、最近の研究開発分野で用いられている「密度汎関数法」のことまで、幅広く紹介されているバランスの良い教科書
計算化学をもっと深く学ぶ
分子力場を用いる分子動力学シミュレーションから、電子状態計算まで、コンピュータを用いて分子を解析する様々な手法がバランス良く紹介されています
量子化学計算の概要
量子化学計算とは？
量子化学計算は、物質を特徴付ける微視的な性質を知りたいと思うときに、対象とする物質の分子レベルでのふるまいを明らかにすることで、問題を解決するための手がかりを与えてくれる便利な道具です
機能性材料 を設計しようとしたり、新しい薬 を開発しようとするときに、研究・開発の現場では様々な課題に直面することもあります
量子化学計算は、たとえば、物質を構成している分子の形、安定性、分子同士の相互作用の強さ、化学反応の起こりやすさなどについて、ある程度の定量的な精度で信頼できる予測を理論的に導きだすことができます
量子化学計算を使ってできること
分子の立体構造を予測できます
https://gyazo.com/8138f25ffada153a8f7e0a3f28c543ca
分子の振動状態（赤外・ラマンスペクトル）を予測できます
https://gyazo.com/48f5e3d868003070ac0fbf8dd02f1e06
分子の励起状態（吸収・発光スペクトル）を予測できます
https://scrapbox.io/files/64e5b1c4f360e5001ba2f16f.webp
分子の電子状態に基づいて様々な解析ができます
https://gyazo.com/347521b2da6cb5361d0b0bbf900fcd87
分子の化学反応を予測できます
https://gyazo.com/95d117bed9703149d3e06b4a5dd67a10
分子シミュレーション手法の比較
https://gyazo.com/4de55b9d707c79158245d1f908545fb6
各アプリの公式 Web サイト
孤立分子 ➡ 量子化学計算
固体物理
第一原理計算のデータを深層学習したニューラルネットワークポテンシャル（NNP）力場を用いる
生体分子
AmebrTools は無償配布
量子化学計算の代表的なソフトウェア
業界標準の Gaussian
Pople は、いろいろあって、Gaussian の使用を拒否されてしまったようです 名前は、分子軌道を表現するために用いる ガウス関数（Gaussian） に由来
ガウス関数：$ \exp(-\alpha r^2)
スレーター関数：$ \exp(-\gamma r)
GaussView というアプリと併せて使うと、高機能な量子化学計算をお手軽に実行できる
https://i.gyazo.com/8ce4153a5a97629bd004dd4d55ccebca.png
計算化学者だけではなく、実験化学者にも Gaussian ユーザーが多い
分からないこと・困ったことがあるとき、適切なアドバイスをしてくれるユーザーがすぐに見つかる
インターネットの検索で有益な情報をたくさん得ることができる
最新版は Gaussian 16
エネルギー計算 や 構造最適化 の収束性が良い
yamnor.icon 山本は、他の量子化学計算プログラムを使うときにも、構造最適化だけは Gaussian で実行する ことがあります
Gaussian の入力ファイル
例：エチレン（$ \mathsf{C_2H_4}）の構造最適化
code:Gausian
C2H4
0 1
C 0.00000000 0.65855185 0.00000000
H 0.91494130 1.22519288 0.00000000
H -0.91468738 1.22560572 0.00000000
C 0.00000000 -0.65855185 0.00000000
H -0.91494130 -1.22519288 0.00000000
H 0.91468738 -1.22560572 0.00000000
実演：Gaussian/GaussView を用いたエチレンの量子化学計算
https://scrapbox.io/files/63f98ded457132001dc5b0e8.mov
Gaussian を販売している国内代理店
無料で使える GAMESS
GAMESS は、Gaussian に次ぐ規模での利用実績がある
General Atomic and Molecular Electronic Structure System の頭文字
アメリカ版（GAMESS-US）、イギリス版（GAMESS-UK）、Firefly という３つの異なるバージョンがある
GAMESSと言うとき、アメリカ版（GAMESS-US） を指すことが多い（yamnor.icon と思います）
どのバージョンも、1970年代後半にアメリカ・エネルギー省の NRCC（National Resource for Computational Chemistry） というプロジェクトで開発されたプログラムコードが共通の先祖
1981年
オリジナル版 GAMESS が アメリカ版（GAMESS-US） と イギリス版（GAMESS-UK） に枝分かれ
1997年
アメリカ版（GAMESS-US）を元にして ロシアで Firefly の開発がはじまった
この３つの異なるバージョンの GAMESS は、別々のグループで開発されている
各グループで独立に、より効率的に計算できるように改良、新しい機能が追加
現在では、お互いに全く異なる性質のものに
アメリカ版 GAMESS
学術＆商用で 無償利用可
➡ 商用でも無償という点は大きい！
オープンソース ではない
プログラムの再配布やソースコードの再利用は禁止
GAMESS-US を利用した研究成果を公表する場合、開発者らによる下記の 原著論文を引用する G. M. J. Barca, et al., J. Chem. Phys., Vol. 152, 154102 (2020)
アップデートが頻繁なので、バージョン番号ではなく、更新日で管理 されている
最先端で開発が進められている 様々な新しい手法 を実行できる
ソースコード を無償で入手することができる
オリジナルの機能を自分で追加する ことができる
GAMESS-US の入力ファイルの例
エチレン（$ \mathsf{C_2H_4}）の構造最適化
code:GAMESS
$CONTRL
SCFTYP=RHF
RUNTYP=OPTIMIZE
ICHARG=0
MULT=1
COORD=CART
$END
$BASIS
GBASIS=N31 NGAUSS=6 NDFUNC=1
$END
$DATA
C2H4
C1 1
C 6 0.00000000 0.65855185 0.00000000
H 1 0.91494130 1.22519288 0.00000000
H 1 -0.91468738 1.22560572 0.00000000
C 6 0.00000000 -0.65855185 0.00000000
H 1 -0.91494130 -1.22519288 0.00000000
H 1 0.91468738 -1.22560572 0.00000000
$END
この講座では
今回は主に、最も利用者が多い Gaussian について説明します
GAMESS の「商用利用でも無償」というメリットは大きいのですが、Gaussian に比べると、利用者は少ないように思います
GAMESS に言及する場合には、３種類のなかで最も代表的なアメリカ版（GAMESS-US）について説明します
量子化学計算を実行できるスパコン
学術利用のみです
産業利用枠もあります
産業利用をメインにした計算環境です
実演：スパコンを用いたエタノールの量子化学計算
https://scrapbox.io/files/63fd3bf4fe79f6001cc36855.mov
量子化学計算を体験できる Web アプリの MolCalc
https://i.gyazo.com/6f67aebc07d389094079c0bbdec6c94a.png
量子化学エンジンは GAMESS
計算レベルは HF/STO-3G
計算手法：Hartree-Fock 法
基底関数：STO-3G
分子を SMILES 記法で指定できる
★ SMILES 記法
SMILES記法とは？
分子の化学構造を 文字列化 して表現する表記方法
Simplified Molecular Input Line Entry System の頭文字
基本的なルール
原子は 元素記号 で表し、水素原子 は 省略 する
二重結合 は「=」、三重結合 は「#」で表す
単結合 は（通常は）省略 する
環構造 は、つながっている原子の後ろに数字でラベル付け する（C1 など）
例：プロパン：CCC、シクロプロパン：C1CC1
芳香環 を構成する原子は小文字にする（c など）
例：ベンゼン：c1ccccc1
table:記述例１
化合物名 化学式 SMILES
メタン CH4 C
アンモニア NH3 N
水 H2O O
二酸化炭素 CO2 O=C=O
窒素 N2 N#N
酸素 O2 O=O
エタン C2H6 CC
エチレン C2H4 C=C
アセチレン C2H2 C#C
table:記述例２
化合物名 化学式 SMILES
ブタン C4H10 CCCC
1,3-ブタジエン C4H6 C=CC=C
1,2-ブタジエン C4H6 C=C=CC
シクロブタン C4H8 C1CCC1
シクロブタジエン C4H4 C1=CC=C1
シクロヘキサン C6H12 C1CCCCC1
シクロヘキセン C6H10 C1CCC=CC1
1,4-シクロヘキサジエン C6H8 C1C=CCC=C1
ベンゼン C6H6 c1ccccc1
実演：MolCalc を用いたエチレンの量子化学計算
https://scrapbox.io/files/63f98604ec8532001c7dcf9c.mov
SMILES 記法で分子を指定する
C=Cと入力する
分子が表示されたら、Calculate Properties をクリック
Indeed をクリック
量子化学計算がはじまる
計算が終わったら、調べたい性質のタブをクリック
Free Energies ：熱力学物性
Vibrational Frequencies：振動解析
Molecular Orbitals：分子軌道
Polarity and Solvations：静電ポテンシャルや双極子モーメントなど
Gaussian や GAMESS を便利に使うための WebMO
Web アプリ
ブラウザーから量子化学計算を実行できる
https://i.gyazo.com/daed19cf1c4a8244ff74da742678b451.png
様々な量子化学計算パッケージのインターフェースとして使える
Gaussian / GAMESS / Molpro / Mopac / NWChem / ORCA / PSI4 / Quantum Espresso / Q-Chem / VASP
Username: guest
Password: guest
30秒までの計算が可能
タイミングが悪いと混んでいる
分子を SMILES で指定するときには、メニューから...
Lookup → Import by → SMILES
実演：WebMO を用いたエチレンの量子化学計算
https://scrapbox.io/files/63f99217e5f976001c763dc8.mov
量子化学計算のコスト
量子化学計算プログラム
Gaussian
GAMESS
学術 / 商用：0 円
GUI アプリ
GaussView
学術：約 16 万円
商用：約 30 万円
学術：1,200 USD（約 16 万円）
商用：3,000 USD（約 41 万円）
計算用サーバー
約 17 万円 / 年
量子化学計算の代表的な計算方法
代表的な計算方法
電子相関を考慮しない計算方法
Hartree-Fock（HF）法
電子相関を考慮する計算方法
密度汎関数法
★ Møller-Plesset 摂動法（例：MP2, MP4）
★ 結合クラスター法（例：CCSD, CCSD(T)）
★ 配置間相互作用（例：CISD, Full CI）法
代表的な計算方法の精度とコスト
https://gyazo.com/f0ec4c466ed23d1325f6ecef79d54eac
Hartree-Fock 法とは？
シュレディンガー方程式
$ \left[ -\frac{1}{2} \sum_i^n \nabla_i^2 + \sum_i^n v(\mathbf{r}_i) + \sum_{i}^n \sum_{j > i}^n \frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|} \right] \Psi = E \Psi
２校目の$ v(\mathbf{r}_i)は、$ i番目の電子に外部（原子核）からはたらくポテンシャル
$ v(\mathbf{r}_i) = -\sum_a^{N_\mathbf{nuclei}} \frac{Z_a}{|\mathbf{R}_a - \mathbf{r}_i|}
３項目の $ \frac{1}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}は、$ i番目と$ j番目の「電子間にはたらく相互作用（電子相関）」を表す
「ある電子は、他の電子の位置に依存しながら運動している」という描像（多電子描像）
➡ 厳密に取り扱うことが難しいため、何らかの「近似」が必要
https://gyazo.com/1ccadafe448870f18efb9fc919e8b028
Hartree-Fock (HF) 方程式
「ある電子は、他の電子の作る平均化された場の中で（他の電子と）独立に運動している」という描像を考える
➡「一電子近似」と呼ぶ
量子化学計算の出発点となる最も基本的な近似
具体的には、$ n電子系の波動関数 $ \Psi(\xi_1,\xi_2,...,\xi_n)を変数分離して、適切な関数$ \psi(\xi_i)の単純な積（Hartree 積）で表すことで、それぞれの電子を１電子軌道$ \psi(\xi_i)に収容する
$ \Psi(\xi_1,\xi_2,...,\xi_n) = \psi_1(\xi_1)\psi_2(\xi_2)\cdots\psi_n(\xi_n)
➡ この「軌道の単純な積」は、Pauli の原理（任意の２つの電子の座標$ \xi_iと$ \xi_jの交換に対して反対称関数であること）を満たさないので、実際には、Slater 行列式で表す必要がある
１電子軌道を用いると...
$ \left[ -\frac{1}{2} \nabla_i^2 + v(\mathbf{r}_i) + V_i^{\mathrm HF} \right] \psi_i = \epsilon_i \psi_i
３項目の$ V_i^{\mathrm HF}は、$ i番目の電子と「他の電子が作る「平均場」に由来する反発エネルギー」
Hartree-Fock 法の計算例（He 原子の場合）
ステップ①：対象とする分子系のハミルトニアン（核座標の配置 $ \mathbf{R}_\alpha、核電荷 $ Z_\alpha、電子数 $ nなどの情報）、初期の参照軌道 $ \Psi^{0}を設定する
https://gyazo.com/b58e2ab39eb27ec51b2959c5cccb8772
対象とする分子系のハミルトニアン：
$ \hat{H} = -\frac{1}{2} \nabla_1^2 -\frac{1}{2} \nabla_2^2 - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{r_{12}}
初期の参照軌道：
$ \Psi^{0}(\xi_1,\xi_2)=\frac{1}{\sqrt{2}}\left[ \alpha(s_1)\beta(s_2) - \alpha(s_2)\beta(s_1)\right]\psi^{0}(\mathbf{r}_1)\;\psi^{0}(\mathbf{r}_2)
ステップ②：「電子２がつくる平均場」中にある電子１について、HF 方程式を解く
https://gyazo.com/646d1e352573e062878bc1e6c273adb7
$ \left[-\frac{1}{2}\nabla_1^2 - \frac{2}{r_1} + \int d\mathbf{r}_2\;|\psi^0(\mathbf{r}_2)|^2 \,\frac{1}{r_{12}}\right]\psi(\mathbf{r}_1)=\epsilon_1 \psi(\mathbf{r}_1)
ステップ③：「電子１がつくる平均場」中にある電子２について、HF 方程式を解く
https://gyazo.com/97545e45b0850511b8931990866f69fa
$ \left[-\frac{1}{2}\nabla_2^2 - \frac{2}{r_2} + \int d\mathbf{r}_1\;|\psi^0(\mathbf{r}_1)|^2 \,\frac{1}{r_{12}}\right]\psi(\mathbf{r}_2)=\epsilon_2 \psi(\mathbf{r}_2)
ステップ④：ステップ②と③で得られた軌道エネルギー $ \epsilon_iと１電子軌道 $ \psi(\mathbf{r}_i)から、全電子エネルギー$ Eを求める
得られた全エネルギー $ E_\mathrm{new}と前回のもの $ E_\mathrm{old}とのエネルギー差 $ \Delta Eを比べたときに...
$ \Delta Eが十分に小さい
➡ 計算終了
$ \Delta Eが大きい
➡ アップデートした分子軌道 $ \Psi_\mathrm{new}を参照軌道 $ \Psi^0に使って、ステップ②〜④を再度実行
https://gyazo.com/9a1f315a77516ed00810e2689f098eb2 https://gyazo.com/42983abf5101a0c917fd57a8443440a2
Hartree-Fock 法の特徴
分子の立体構造や化学的性質について、多くの場合、十分な精度で予測できる
分子が最適構造にあるとき、その分子の全エネルギーの値を99%程度の精度で見積もることができる
しかし、化学反応を議論するためには、99%程度という精度でも不十分であり、より高精度な予測が必要となる
化学反応に伴う活性化エネルギー $ E_aを議論しようとする場合、$ E_aは 20 kJ/mol 程度（5 kcal/mol 程度）になることが多いので、この程度以下の誤差となることが望ましい
電子相関を考慮するように HF 法を拡張し、HF 法以上の計算精度を持つ手法のことを「post Hartree-Fock 法」と呼ぶ
post HF 法では、HF 法では無視してしまった「電子相関」の効果について、様々な工夫で考慮する
電子相関とは？
Hartree-Fock 法は、電子同士の相互作用（電子相関）を平均化して扱うという近似（平均場近似）をしている
https://gyazo.com/1ccadafe448870f18efb9fc919e8b028
電子相関の大きさは「電子相関エネルギー」（$ \Delta E_\mathrm{corr}）と呼ばれ、一般に、シュレディンガー方程式の厳密解（$ E_\mathrm{exact}）と HF 計算で得られる値（$ E_\mathrm{HF}）との差分として定義される
$ \Delta E_\mathrm{corr} = E_\mathrm{exact} - E_\mathrm{HF}
電子相関の具体例
水分子１個の場合で、動的電子相関エネルギーは 140 kcal/mol 程度
動的電子相関を考慮する計算方法を用いると、動的電子相関に由来する誤差は...
MP2 摂動法では 8 kcal/mol 程度
CCSD(T) レベルでは 1 kcal/mol 以下
化学反応を定量的に議論することが可能な精度
電子相関を考慮する計算方法
密度汎関数法
★ Møller-Plesset 摂動法
二次の Møller-Plesset (MP2：second-order Møller-Plesset) 摂動法は、計算コストを比較的抑えながら電子相関の効果を大雑把ながらも比較的精度良く見積もることができる
生体分子など、大規模な分子系に対する電子相関補正として用いられることが多い
★ 結合クラスター法
化学反応にともなうエネルギー変化に基づいて反応速度を議論したい場合など、高い定量性が求められる応用研究に多く利用されてるのは、結合クラスター法を用いた CCSD(T) レベル
ただし、計算手法が高度になるほど、計算するためのコスト（計算に必要な時間やコンピュータの性能）も大きくなってしまうことがほとんどなので、目的に見合った精度とコストのバランスを考えることも大事になる
密度汎関数法（DFT：Density Functional Theory）とは？
近年、研究開発の最前線で大活躍している計算手法
電子密度を試行関数とする多電子状態問題の近似解法
https://gyazo.com/f3b401fbb48d795614ae3ef5ac4b808c
電子相関の寄与について、 Hartree-Fock 法とほぼ同程度の計算コストで考慮することができる
Hartree-Fock 法 などは波動関数を試行関数とする方法であるが、DFT は電子密度を試行関数とする方法であり、その基盤となる考え方が大きく異なっている
HF 法の基礎は HF 方程式だが、DFTを実践する上で鍵になるのはコーン・シャム（Kohn-Sham; KS）方程式
Kohn-Sham 方程式とは？
HF 方程式
$ \left[ -\frac{1}{2} \nabla_i^2 + v(\mathbf{r}_i) + V_i^\mathrm{HF} \right] \psi_i = \epsilon_i \psi_i
３項目の$ V^\mathrm{HF}は「ある電子-他の電子間の「平均」反発エネルギー」を表す
「ある電子は、他の電子の作る平均化された場の中で（他の電子と）独立に運動している」という描像
KS 方程式
$ \left[ -\frac{1}{2} \nabla_i^2 + v(\mathbf{r}_i) + V_i^\mathrm{eff} \right] \psi_i = \epsilon_i \psi_i
３項目の $ V^\mathrm{eff}は、「外場有効ポテンシャル」
$ V^\mathrm{eff} = J + V_\mathrm{xc}
$ J：電子-電子間の古典的な静電（クーロン）相互作用
$ V_\mathrm{xc}：非古典的な相互作用であり、（交換相関）汎関数と呼ばれる
➡ DFT を特徴付ける重要なもの
KS 方程式とHF 方程式の違いは、３項目が「平均場$ V^\mathrm{HF}」か「外場有効ポテンシャル$ V^\mathrm{eff}」かのみ
つまり、KS 方程式は、HF 方程式と同程度の時間で計算できる
HF 方程式中の「平均場$ V^\mathrm{HF}」は「既知の関数」であるが、KS 方程式中の「交換相関汎関数$ V_\mathrm{xc}」は厳密な形が不明な「未知の関数」である
$ V_\mathrm{xc}の厳密な形が分からないので、人為的に仮定された汎関数を使う
この「人為的に仮定された」汎関数の質が、DFT の精度を決定することになる
密度汎関数法で用いる「汎関数」
汎関数とは？
DFT の計算精度と信頼性は、電子密度と系のエネルギーを関係付ける汎関数の種類に依存する
適切な汎関数を選べば、良い計算精度が得られる
交換・相関汎関数の分類
局所スピン密度近似（LSDA）型
電子密度 $ \rho のみで表現された汎関数
一般化勾配近似（GGA）型
LSDA を電子密度勾配 $ \nabla \rho で補正した汎関数
meta GGA 型
GGA を運動エネルギー密度 $ \tauで補正した汎関数
hybrid 型
Hartree-Fock 交換積分$ E_\mathrm{x}^\mathrm{HF}を一定量混合した汎関数
https://gyazo.com/f9a0046baad0a2b25cf7d682c8f9adbb
交換・相関汎関数の種類と計算精度
https://gyazo.com/4b97ecc378e9803c33fe1f464224cda2
交換・相関汎関数の階層的な分類
Perdew によって提案された汎関数の階層的な分類です。
https://gyazo.com/70fdb4eec4ffaa4eee9e5d8643beecaf
交換・相関汎関数の例
代表的な200個の交換・相関汎関数をまとめると...
Local：非 hybrid 型
Disp：分散力補正あり
下線：長距離補正あり
https://gyazo.com/7f3a9a074aae8c94477556d2bb3c3da5
いくつかの代表的な汎関数
基本的な GGA 汎関数
BLYP （Becke-Lee-Yang-Parr）
PBE（Perdew-Burke-Ernzerhof）
meta GGA 汎関数
BMK（Boese-Martin functional for Kinetics）
ミネソタシリーズ（M06 など）
Hartree-Fock 交換項と組み合わせた hybrid 型の汎関数
B3LYP
以前は実質上の業界標準（デファクトスタンダード）であるかのように多用
長距離電子相関の見積もりが不得意であることが知られるようになった
長距離補正を含む「CAM-B3LYP」、分散力補正を含む「B3LYP-D3」などが開発されている
PBE0
PBEPBE
長距離補正を含む汎関数
CAM-B3LYP
ωB97シリーズ（ωB97, ωB97-X, ωB97-XD）
LC-BOP
分散力補正を含む汎関数
B3LYP-D3
ωB97X-D
ミネソタシリーズ（M06-2X, X11 など）
hybrid 交換・相関汎関数
hybrid 型は、LSDA / GGA / mGGA などの交換汎関数に対して、一定の割合で Hartree-Fock 交換積分$ E_\mathrm{x}^\mathrm{HF} を混合した汎関数です。
非 hybrid 型（local 型）と hybrid 型の汎関数ペア（例：BLYP vs B3LYP）について、計算される物性値の平均二乗偏差（RMSD）を比較してみると...
NCD：分子間相互作用エネルギー
ID：異性化エネルギー
BH：反応障壁の高さ
https://gyazo.com/eecd63c52ff0735df7ece58b720e7b39
$ E_\mathrm{x}^\mathrm{HF}を適度に hybrid することで計算精度が向上
分散力補正
標準的な DFT の弱点のひとつは、分散力（van der Waals 相互作用）の記述ができない
➡ このために、分子集合体における分子間相互作用を過小評価する傾向にある
Grimme は、分子力場などで用いられるように経験的なパラメータを用いて、各原子対に対する引力的なエネルギー項を付加する分散力補正を提案
Grimme の分散力補正のあり・なし（例：BLYP と BLYP-D2, BLYP-D3）について、計算される物性値の平均二乗偏差（RMSD）を比較してみると...
NCED：分子間相互作用エネルギー
IE：異性化エネルギー
https://gyazo.com/260c089bbd967b6b9b70a2d5f55d930e
Grimme の分散力補正を加えることで、コストを追加することなく、計算精度を向上させることができる
交換・相関汎関数のベンチマーク
200 個の汎関数について、ベンチマーク計算をしてみると...
汎関数名の右にある数字は、計算精度のランキング
左にあるラベルは、評価した物性値：
NCED / NCED / NCD：分子間相互作用エネルギー
IE / ID：異性化エネルギー
TCE / TCD：熱化学量
BH：反応障壁の高さ
https://gyazo.com/42c648b2b9ebf15b3674168b8b6de3fc
GGA から meta GGA に変えても、多くの場合、計算精度はそれほど良くはならない
非 hybrid 型（local 型）から hybrid 型に変えると、ほとんどの場合、計算精度が向上
量子化学計算の代表的な基底関数
基底関数とは？
任意の波動関数 $ \Psi を有限個の既知の関数 $ \phi_1, \phi_2, ..., \phi_n の線形結合で表すという近似（線形近似）を考える
$ \Psi = c_1 \phi_1 + c_2 \phi_2 + ... + c_n \phi_n = \sum_i^n c_i \phi_i
この既知の関数 $ \phi_1, \phi_2, ..., \phi_n のことを 基底関数 と呼ぶ
展開係数 $ c_1, c_2, ..., c_nの値を変えると、波動関数 $ \Psiが変わる
量子化学計算では、あらかじめ適切にパラメータを設定した基底関数 $ \{\phi_n\}を使って、展開係数 $ \{c_n\}の値について、最適なものを探そうとする
つまり、計算の途中で...
基底関数 $ \{\phi_n\}は変えずに、
展開係数 $ \{c_n\}の値だけを変える。
https://i.gyazo.com/a4e3310a0265db15ba8d23dbb445192c.png
基底関数には様々な種類があり、どのような分子を解析したいのか、何を調べたいと思っているのか、どのような計算方法を用いるのか、それぞれに応じて、適切に選択する必要がある
基底関数の個数が多いほど、電子分布を柔軟に表現できるので、計算精度が高くなる
計算精度が高くなればなるほど、計算コストも高くなってしまう
計算精度と計算コストのバランスを考えて、適切に選ぶ
水素分子の場合
水素分子を表す波動関数 $ \Psi は、たとえば、２個の水素原子 $ \mathrm{H}_aと$ \mathrm{H}_bの中心に置いた1s 原子軌道$ \phi_aと$ \phi_bの線形結合として、次のように表すことができる
$ \Psi = c_a \, \phi_a + c_b \, \phi_b
水素原子については、1s 軌道の厳密解が次のように分かっている
$ \psi_\mathrm{H1s} = \sqrt{\frac{1}{\pi}} \exp(-r)
https://i.gyazo.com/dcfa9f3726de4f8d2c15e504cbece8b1.png
したがって、水素原子の 1s 軌道を参考にして、$ \exp(-\zeta r)型の関数を基底関数として使っても良さそうだ
$ \phi_a = \exp(-\zeta r)
このような関数を スレーター型関数 と呼ぶ
係数 $ \zetaは、関数の「広がりの程度」を表す変数
水素原子の厳密な解析解（= 電子の分布を適切に表すことができる）なので、計算精度は高い
https://gyazo.com/0857066e42c0815178e589a5802444d1
https://gyazo.com/dbb67c916e73c77f44cf84a176b59e55
ただし、数値計算の計算効率は悪い
yamnor.icon 図にあるように、原子の中心（$ r = 0）で値が不連続になるから、数値的に積分することが難しい
計算効率が悪いので、ふつうは、量子化学計算では使わない
ガウス型関数
量子化学計算では、スレーター型関数$ \exp(-\zeta r)の代わりに、$ \exp(-\alpha r^2)型の関数を基底関数として使う
$ \phi_a = \exp(-\alpha r^2)
このような関数を ガウス型関数 と呼ぶ
https://gyazo.com/66a23423bf1cc335cf3e817cdf51f374
係数 $ \alphaは、関数の「広がりの程度」を表す変数
量子化学計算の前に、あらかじめ設定されている定数（計算の途中で変化させない）
ガウス型関数は、原子軌道に対する近似解なので、計算精度は高くはない
ただし、数値計算の計算効率が良い
yamnor.icon 上の図にあるように、原子の中心（$ r = 0）でも値が連続になるから、数値積分が簡単
計算精度の悪さは、拡がりの程度を変えたガウス型の基底関数をいくつか線形結合させることでカバーする
いくつかのガウス型関数を線形結合させたものを短縮型ガウス基底とよぶ
たとえば...
$ \phi_\mathrm{赤} = c_\mathrm{緑} \, \chi_\mathrm{緑} + c_\mathrm{紫}\, \chi_\mathrm{紫} + c_\mathrm{黄}\, \chi_\mathrm{黄}
$ = 0.5 e^{-0.5 r^2} + 0.3 e^{-1.0 r^2} + 0.2 e^{-0.3 r^2}
https://gyazo.com/40534f7c97701c1c01c4e82ecef357c2
赤い実線で示す短縮型ガウス関数は、青い点線で示すスレーター型関数のかたちを、まあまあ良く再現している
基底関数の種類①
最小基底系（STO-nG基底系）
例
STO-3G
３種類のガウス型関数（3G）を足し合わせて、スレーター型軌道（STO）に似せたもの
STO-4G
STO-5G
計算精度の比較
STO-3G < STO-4G < STO-5G
計算コストは低いが、計算精度も低いので、実験結果と比較するなどの実用的な解析に用いることはあまりない
基底関数の種類②
Pople 系
例
3-21G
6-31G
6-311G
計算精度の比較
3-21G < 6-31G < 6-311G
最小基底系に比べると、より柔軟に電子のふるまいを記述できるため、計算精度が改善される
原子価軌道に対して、拡がりの異なる複数の関数を組み合わせて表現する
一般的な有機化合物に対しては、基盤となる 6-31G に 分極関数 を加えた 6-31G(d) や 6-31G(d,p) などが、計算効率が良い・実用的な精度をもつ基底関数としてよく用いられている
分割価電子（split valence）基底関数
内殻軌道には１つ、価電子軌道には複数の短縮ガウス型関数を用いた基底関数系
Pople 系の場合
https://gyazo.com/262c6f95b97dc36484421eb601dbb708
6-31Gの場合
１つの占有軌道に対して...
内殻軌道
６個のガウス型を足し合わせた（短縮した）１種類の関数で記述
価電子軌道
３個のガウス型を短縮したもの＋別の１個のガウス型を組み合わせて、計２種類の関数で記述
２種類の関数で原子価軌道を表すものを double zeta (DZ) 型とよぶ
6-311Gの場合
１つの占有軌道に対して...
内殻軌道
６個のガウス型を足し合わせた（短縮した）１種類の関数で記述
価電子軌道
３個のガウス型を短縮したもの＋別の１個のガウス型＋別の１個のガウス型を組み合わせて、計３種類の関数で記述
３種類の関数で価電子軌道を表すものを triple zeta (TZ) 型とよぶ
具体例
https://gyazo.com/5dd401bc84a174ea0e3e2aff0b703f35
基底関数の種類③
高精度計算に適した基底関数（Correlation-Consistent 基底系）
計算精度の比較
cc-pVDZ < cc-pVTZ < cc-pVQZ
応用研究の中では、double zeta 精度の cc-pVDZ が頻繁に用いられる
double zeta 精度では、電子相関効果の取り込みが不十分である場合もある
化学反応などにおいて 10 kcal/mol 程度のエネルギー差を定量的に議論したいときには、cc-pVTZ 程度（triple zeta 精度）の基底関数を CCSD(T) レベルの計算手法と組み合わせて用いる
分極関数とは？
分子の複雑な構造を適切に記述したい場合には、各原子上に配置した基底関数に 分極関数（polarization function） を追加することで、電子分布が柔軟に変形する自由度を与えることができる
たとえば...
p 軌道をもつ水素以外の原子に対しては、d 型の関数を追加する
http://4.bp.blogspot.com/-9odFiZctEe0/VX1_xMr2R6I/AAAAAAAADi4/uLpgZkL2MU0/s1600/Fig2.7.8.jpg
球対称な s 軌道のみを持つ水素原子に対しては、p 型の関数を追加する
http://2.bp.blogspot.com/-nQJzi068Llo/VX1_kBl-cZI/AAAAAAAADiw/aE4Ln8gsYFM/s320/Fig2.7.7.jpg
表記方法 / 指定方法
Pople 系に対して、分子内にある水素以外の原子に分極関数を追加する場合
3-21G(d) or 3-21G*
6-31G(d) or 6-31G*
Pople 系に対して、水素以外の原子に分極関数を追加した上で、水素原子にも分極関数を追加する場合
3-21G(d,p) or 3-21G**
6-31G(d,p) or 6-31G**
cc-pVnZ 基底関数（$ n = D, T, Q, 5, 6）は、水素原子および水素以外の原子のどちらにも、全ての原子に分極関数が追加される
分散関数とは？
負電荷を帯びた分子は、中性分子と比較して、電子分布が広がっている場合がある
このような場合、物性諸量を正しく見積もるためには、電子分布の広がりを記述するために 分散関数（diffuse function） を加えることができる
表記方法 / 指定方法
Pople 系に対して、分子内にある水素以外の原子に分散関数を追加する]場合
3-21+G(d) or 3-21+G*
6-31+G(d) or 6-31+G*
Pople 系に対して、水素以外の原子に分散関数を追加した上で、水素原子にも分散関数を追加する場合
3-21++G(d,p)
6-31++G(d,p)
cc-pVnZ 基底関数（$ n = D, T, Q, 5, 6）に対して、水素原子および水素以外の原子のどちらにも、全ての原子に分散関数を追加する場合
aug-cc-pVDZ
aug-cc-pVTZ
基底関数のベンチマーク（VSCF-PT2 計算の場合）
https://gyazo.com/bc41c8ddebc37071865e26be60e1490a
基底関数と計算時間（VSCF-PT2 計算の場合）
https://gyazo.com/b19c75dc1b69aaddb5d92f24829e676e
オススメの計算手法と基底関数は？
yamnor.icon は、一般的な有機化合物の基底状態の物性を調べようとする場合、まずは B3LYP/6-31G(d) や B3LYP/6-31G(d,p) を使うことが多いです
原子間の結合距離を定量的に議論するには、最低でも double zeta 精度（6-31G）が必要
原子間の結合角度を定量的に議論するには、分極関数（6-31G(d)）が必要
DFT 汎関数と基底関数の組み合わせ
https://gyazo.com/bfb3f8fdd4fe483d734df2659bd6a553
計算条件の表し方
量子化学計算の計算精度や計算コストは、主に、計算方法と基底関数の組み合わせで決まる
計算方法の種類 と 基底関数の種類 を スラッシュ（/） で挟んで表記することが多い
たとえば...
B3LYP/6-31G(d)
MP2/cc-pVDZ
あるレベルで 構造最適化 した後、より高精度なレベルで エネルギー計算 をおこなった場合、２つの計算レベルを ダブルスラッシュ（//） で挟んで表記することがある
たとえば...
CCSD(T)/aug-cc-pVTZ//B3LYP/6-31G(d)
最初に、手軽に電子相関の寄与を考慮できる密度汎関数法（B3LYP汎関数）と小さな基底関数（6-31G(d)）の組み合わせで、構造最適化を実行
次に、B3LYP/6-31G(d)レベルで最適化した分子構造を用いて、高精度な量子化学計算手法（CCSD(T)）と大きな基底関数（aug-cc-pVTZ）の組み合わせで、エネルギー計算を実行
量子化学計算に基づく構造最適化
分子の安定構造の探索
各原子核に働く力（エネルギー勾配 $ \frac{\partial E}{\partial Q}）の方向にカタチをちょっとずつ変化（$ x_\mathrm{new} = x_\mathrm{old} - \alpha \frac{\partial E}{\partial Q}）させることで、分子はより安定となる（ = エネルギーが小さくなる）
$ E(x_\mathrm{new}) < E(x_\mathrm{old})
https://gyazo.com/6fad80947159d295366c18cfb474210c
「カタチをちょっとずつ変化させる」を繰り返すことで、分子の立体構造を半自動的に最適化できる
実際には、多くの量子化学計算プログラムは、エネルギー勾配の方向に動かすという単純な方法（勾配降下法）ではなく、もっと効率的なアルゴリズムを使っている
Gaussian では、デフォルトでは、GEDIISというアルゴリズムを用いる
実演：Gaussian を用いた水分子の安定構造の探索
https://scrapbox.io/files/63f9a0bcd06e52001c9d4a71.mov
最適化される様子を確かめるため、HOH 結合角を 150° 程度にしたものを初期構造とした
計算例：水分子の構造最適化
https://gyazo.com/c2e381a0aff6f0bfb93bac93b78c270f
量子化学計算に基づく振動解析
振動解析とは？
分子固有の振動数と赤外強度を理論的に計算する方法
実験的に測定された赤外スペクトルの解釈を補助する手段として活用されている
各振動モードに対するラマン強度を計算することも可能
赤外吸収スペクトルとは？
水分子の赤外吸収スペクトル（実験結果）
https://gyazo.com/246c01937cec7d08b6c54b43d41e6f3e
分子振動とは？
水（$ \mathrm{H_2O}）の変角振動 / 対称伸縮振動 / 逆対称伸縮振動
https://gyazo.com/e4d9bb89f0f61f0533b9532d0e024287https://gyazo.com/145fab5f0795498d91045fd285f6b556https://gyazo.com/80f8eb416e3294f461be22f863b7cc58
振動解析の計算手順
対象とする分子について、構造最適化を実行（安定構造を探索）する
構造最適化を実行したときと同じ量子化学計算レベルで、振動解析を実行する
実演：Gaussian を用いた水分子の振動解析
https://scrapbox.io/files/63f9de9406a787001b9cc00a.mov
計算例：水分子の振動解析
https://gyazo.com/288677666289ced8a6fec15061e4bcb8
分子振動の非調和性
高い精度（MP2/cc-pVTZ）で計算した場合
https://gyazo.com/e1d212dc0d548ac021fe5ffb7701964a
MP2摂動法を用いることで、実験値にかなり近くなる
高波数の振動モード（対称伸縮＆逆対称伸縮）では、誤差が大きい
分子振動を調和振動子だと見なす近似（調和近似）による誤差（非調和性）の影響だと考えられる https://gyazo.com/dcff54c2ebf09e24e2db0b7bfa25ff73
非調和性の問題を手軽に解決する場合には、振動数に適当な数字（スケーリング因子）を掛けて補正する
事例：
HF/6-31G(d) の場合 → 0.903
B3LYP/6-31G(d,p) の場合 → 0.961
https://i.gyazo.com/ec1b28045886c70d80aa644f4f88a593.png
★ 非調和振動数解析
分子振動をより定量的精度で解析したい場合には、調和近似を越えた取り扱いが必要
GAMESS では、非調和性を考慮した振動解析として、Vibrational SCF (VSCF) 法を利用することができる
VSCF法で計算した場合
https://gyazo.com/d66cdc593c263cd07ac4f263e4ba3715
量子化学計算に基づく化学反応の解析
分子の遷移状態
ある化学反応について、反応速度を議論するときには、その反応に伴う遷移状態を調べる必要がある
たとえば、ある反応の速度定数 $ kについて、始状態と遷移状態とのエネルギー差、つまり、活性化エネルギー $ \Delta G^{\ddagger}が分かっていれば、Eyring-Polanyi 式から...
$ k = \chi \frac{k_\mathrm{B}T}{h}\exp\left(-\frac{\Delta G^{\ddagger}}{RT}\right)
$ \chi：遷移状態において、反応座標に沿った分子振動が一往復したときに反応する確率（透過係数）
例：Diels-Alder 反応
https://gyazo.com/6245b572e337bc0273b0b260af89a429
https://i.gyazo.com/a53e519a0a2b7c54971c9b3bd1ef1871.png
遷移状態の探索
多くの量子化学計算プログラムには、遷移状態を探索する機能が備わっている
単純に「構造最適化の逆」（$ x_\mathrm{new} = x_\mathrm{old} + \alpha \frac{\partial E}{\partial Q}）をやっても、うまくいかない
https://gyazo.com/6fad80947159d295366c18cfb474210c
エネルギーの二次微分（$ \frac{\partial^2 E}{\partial Q^2}）を要素にもつ Hessian 行列 の固有値（力の定数 or 曲率：$ \chi）を求めて...
Hessian 行列の負の固有値を持つ固有ベクトルの方向に少しずつ動かすという方法がある
→ Eigenvector Following 法
調和近似（$ E = \kappa Q^2）からの歪みが大きな方向に少しずつ動かすという方法がある
Hessian 行列を計算するためのコストがまあまあ大きい
➡ 構造最適化よりも計算コスト（時間）がかかる
実演：Gaussian を用いた水分子の遷移状態の探索
https://scrapbox.io/files/63f9c6f1a1e4d5001b464477.mov
遷移状態に近い（と予想される）構造から出発する
今回の場合、H-O-H 結合角を大きくする
yamnor.icon 見つけようとしている遷移状態は直線型だと考えられるから
遷移状態探索の設定方法
Job Type
https://i.gyazo.com/c1ba79cb670331d2f93de606d041500f.png
Optimizationを選択
Optimize to a TS (Berry)
Force Constants Calculate at First Point
遷移状態で振動解析すると...
値が負（虚数）となる振動数がある
https://i.gyazo.com/862a1da82d62e9ec3de7e7916cd77573.png
安定構造では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
局所安定状態（下に凸の放物線）であり、エネルギーの座標に関する二次微分の値が常に正である
$ \frac{\partial^2 E}{\partial Q^2} > 0
振動解析すると、すべての振動モードの振動数が正の値になる
https://gyazo.com/e949b48ac2e562a22b8de1963cf0c3ec
遷移状態では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
特定の変位方向に対して局所不安定状態（上に凸の放物線）であり、その変位方向にはエネルギーの座標に関する二次微分の値が負となる
$ \frac{\partial^2 E}{\partial Q^2} < 0
振動解析すると、１つの振動モードの振動数のみが負の値になる
https://gyazo.com/1f0cfca95edab362017ddd028edbef68
実演：Gaussian を用いたブタジエンの cis/trans 異性化反応の解析
https://scrapbox.io/files/63faad6c14c7b1001bf5afac.mov
量子化学計算に基づく電子励起状態の解析
光の吸収とは？
光のエネルギーが分子を基底状態から励起状態に変えるエネルギーとして消費されるプロセス
https://gyazo.com/0b9a593af6e116be42aa973d6e203c8b
$ E_{励起} - E_{基底} = \Delta E = h c / \lambda_{abs}に相当する光のエネルギーを分子が吸収し、基底 → 励起状態への電子遷移が起こる
$ c：光の速度
分子の吸光・発光と基底・励起状態
https://gyazo.com/cc8fd20bf491678066145f87d47dc9a6
吸光
基底（$ \mathrm{S_0}）状態の最適構造 → 励起（$ \mathrm{S_n}）状態
つまり、量子化学計算に取り組むときに、分子の立体構造は、基底状態で最適化する必要がある
発光
励起（$ \mathrm{S_n}）状態の最適構造 → 基底（$ \mathrm{S_0}）状態
つまり、量子化学計算に取り組むときに、分子の立体構造は、励起状態で最適化する必要がある
励起状態の計算方法
TD-DFT 法
時間依存密度汎関数理論（Time-Dependent DFT）法。励起状態の計算方法として、最も普及している
一電子励起配置間相互作用法
Symmetry Adapted Cluster/Configuration Interaction (SAC-CI) 法
完全活性空間 (Complete Active Space) SCF 法
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}
線形応答理論に基づいて、時間依存 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の値が、励起 ＆ 脱励起のエネルギーに対応します。
行列 $ \mathbf{B}を無視して上記の方程式を解く扱いのことを、Tamm-Dancoff 近似（TDA）といいます：
$ \mathbf{AX} = \omega \mathbf{X}
➡ 遷移モーメントの定量性が悪くなる可能性があります。
➡ 一方で、励起エネルギーの計算精度が改善されることがあります。
TD-DFT のベンチマーク（汎関数依存性）
https://gyazo.com/53e8a78345710b4d969c72d291f8d929
領域分割 hybrid 型（長距離補正）を用いることで、計算精度が改善しています。
蛍光（fluorescence）の結果に関しては、ベンチマークに用いたデータが少ないので、信頼性に欠けるかもとのこと。
実演：ホルムアルデヒドの電子励起状態の解析
https://scrapbox.io/files/63faba13bd096f001c1de659.mov
質疑応答
ご質問
さいごに
共同研究・技術相談
計算化学（分子動力学・量子化学・機械学習）に関して、共同研究や技術相談をお受けしています
興味・関心を持っていただけましたら、山本までお問い合わせください
連絡先
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください