oqc2023
https://gyazo.com/7ddbb56165e8abc564c5f750146432fb
資料へのリンク
https://gyazo.com/c3975b18b674ad78e12d91bc1475463e
概要
日時：毎週火曜日１・２限（9 時〜）
場所：653 講義室
はじめに
この資料について
このページでは、山本典史 yamnor.icon が担当する 千葉工業大学 / 大学院 / 工学研究科 / 応用化学専攻 の講義 「有機量子化学特論」 の資料を公開しています
この資料の利用条件については、/about にある注意事項をご確認ください プレゼンテーションモード
資料を共有する時には、プレゼンテーションモードという機能を使っています
ページの右側にある を押して、Start presentationを選ぶと、プレゼンテーションモードになります
23-09-19 第１回
今回の内容
ガイダンス
量子化学計算の概要
ガイダンス
この授業について
この授業では、量子化学計算を使いはじめるときに必要となる知識と技術、量子化学計算とは何か、そしてそれを使うことで何ができるのか、その基本的な概念を説明します。
量子論や量子力学を学んでいない受講生でも、量子化学計算を実行する上で基本となる知識を身につけることができるように、分かりやすく解説します。
具体的な量子化学計算ツールとして、主に GAMESS という、無償で利用できるアプリを中心に紹介します。さらに、可視化ツールである WebMO を活用する方法も詳しく解説します。
受講生がすぐに量子化学計算に取り組めるよう、実際の入力ファイルの作成方法や計算の実行方法、結果の読み取り方なども丁寧に説明します。
この授業では具体的な例を通じて、量子化学計算の応用についても学びます。分子の安定構造、振動状態、励起状態、化学反応の解析など、実際の研究や開発で使われるような問題に取り組みながら、計算手法の実践的な使い方を身につけます。
この授業を受講する受講生には、量子化学計算をすぐに使い始めることができる WebMO のアカウントを発行します。量子化学計算のアプリをインストールする手間なく、量子化学計算を体験できます。
この授業の成績評価
中間試験・期末試験は実施しない
講義内で毎回実施する演習課題（レポート・小テストなど）に基づいて評価をおこなう
点数配分：演習課題（１３回）100％
この授業の参考書
この授業にピッタリの本
授業では GAMESS を使って説明するので、Gaussian を使う場合には、この本も参考にしてみてください
はじめて量子化学計算に取り組む人向けとして、読みやすい入門書
量子化学計算をもっと深く学ぶには
分子力場を用いる分子動力学シミュレーションから、電子状態計算まで、コンピュータを用いて分子を解析する様々な手法がバランス良く紹介されています
量子化学計算の土台となる Hartree-Fock 法について、しっかりと学ぶことができます
量子化学計算でよく用いられている密度汎関数法について、しっかりと学ぶことができます
この授業のスケジュール
ガイダンス + 第１〜９回
担当：山本 典史 yamnor.icon
形式：対面（653 教室）
第１０〜１３回
担当：原口 亮介 先生
形式：オンデマンド（manaba）
量子化学計算の概要
量子化学計算とは？
機能性材料 を設計しようとしたり、新しい薬 を開発しようとするときに、研究・開発の現場では様々な課題に直面することもあります
そんなとき、量子化学計算を活用すると、分子の形、安定性、分子同士の相互作用の強さ、化学反応の起こりやすさなど、対象とする物質の 分子レベルでのふるまい を明らかにすることで、問題を解決するための手がかり を導きだすことができます
量子化学計算を使うと何ができる？
分子の立体構造を予測できます
https://gyazo.com/8138f25ffada153a8f7e0a3f28c543ca
分子の振動状態（赤外・ラマンスペクトル）を予測できます
https://gyazo.com/48f5e3d868003070ac0fbf8dd02f1e06
分子の励起状態（吸収・発光スペクトル）を予測できます
https://scrapbox.io/files/64e5b1c4f360e5001ba2f16f.webp
分子の電子状態に基づいて様々な解析ができます
https://gyazo.com/347521b2da6cb5361d0b0bbf900fcd87
分子の化学反応を予測できます
https://gyazo.com/95d117bed9703149d3e06b4a5dd67a10
【演習課題１】量子化学計算では何ができない？
「量子化学計算ではできないと思うこと」を「５個」考えて、下記の QR コードのリンク先にある「投票」のページから投稿してください。
投稿するときには、「５個」をまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
https://gyazo.com/7b2b5cbecf72ea805b6ae0b3eb99d441 https://gyazo.com/ead8c950c87bb4d4f94073c6f3b78c69
分子シミュレーション手法の比較
https://gyazo.com/4de55b9d707c79158245d1f908545fb6
各アプリの公式 Web サイト
孤立分子 ➡ 量子化学計算
固体物理
第一原理計算のデータを深層学習したニューラルネットワークポテンシャル（NNP）力場を用いる
生体分子
AmebrTools は無償配布
【応用事例】金属有機構造体 MOF でゲスト分子がどのように拡散するのか
LAMMPS を用いた分子動力学シミュレーション
https://gyazo.com/c5d457492844fea35ecca6b1152e74c5
井上＠山本研, 日本コンピュータ化学会（2023）
【応用事例】インフルエンザウイルスの薬剤耐性化はどのようにおこるのか？
Amber を用いた分子動力学シミュレーション
https://gyazo.com/28406ac07980e525221bf2a2e861c2fe
Mohini＠山本研 et al, PeerJ (2021)
【演習課題２】量子化学計算の活用方法
次に示すような「新しい薬を開発するケース」について、量子化学計算を用いることが適切だと思われるものをすべて選んでください。
x 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子を合成するための準備として、合成経路の素反応のひとつについて、対応する遷移状態のエネルギーを知りたい。 x 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子について、分子内でどのような電荷分布になっているのか知りたい。 _ 新しい薬を設計するときに、１００万個の化合物データベースの中から、薬物候補化合物としてふさわしい分子を１０００個選びたい。 x 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子の赤外吸収スペクトルについて、それぞれのピークがどのような振動モードに由来するのかを調べたい。 _ 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子について、どのような結晶構造になっているのかを予測したい。 量子化学計算を体験できる 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：静電ポテンシャルや双極子モーメントなど
【演習課題３】MolCalc を用いたエタノールの量子化学計算
MolCalc を用いてエタノール（$ \mathsf{C_2H_5OH}）の量子化学計算をおこない、この分子について、次に示す化学的性質を予測してください。 生成熱
Free Energies → Heat of Formation
-237.88 $ \mathrm{kJ\;mol^{-1}}
O-H 伸縮モードの振動数
Vibrational Frequencies
3910.9 $ \mathrm{cm}^{-1}
HOMO と LUMO のエネルギー差
Molecular Orbitals
25.05 eV
双極子モーメントの値
Polarity and Solvation → Dipole
1.68 Debye
23-09-26 第２回
今回の内容
量子化学計算の代表的なソフトウェア
量子化学計算の代表的なソフトウェア
業界標準の 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
【演習課題１】量子化学計算の代表的なソフトウェア
量子化学計算の代表的なソフトウェアに関する次の文章について、適切なものをすべて選んでください。
_ 量子化学計算のプログラムは数多くあるが、デファクトスタンダード（事実上の業界標準）だと言えるのは、アイオワ州立大学の John A. Pople らが中心に開発している GAMESS である。 x GAMESS-US は、オープンソースではなく、プログラムの再配布やソースコードの再利用は禁止されているが、学術・商用のどちらの用途でも無償で利用することができる。 x Gaussian は、有料の量子化学計算プログラムであり、GaussView というアプリと併せて使うことで、手軽に実行できる。また、エネルギー計算や構造最適化の収束性が良いことが知られている。 _ GAMESS には、現在、GAMESS-US、GAMESS-UK、Firefly という３つの異なるバージョンがあるが、全く同じ機能が実装されており、計算効率もほぼ同等である。 x ある分子について、Gaussian と GAMESS という異なるプログラムを用いて量子化学計算を実行しても、計算条件が等しければ、全エネルギーの値などについて、ほぼ等しい結果が得られる。 【演習課題２】Gaussian vs GAMESS
量子化学計算の代表的なソフトウェアである「Gaussian」と「GAMESS」の「メリット」と「デメリット」について、下記の QR コードのリンク先にある「投票」のページから投稿してください。
投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
この授業では
今回は主に、最も利用者が多い Gaussian について説明します
GAMESS の「商用利用でも無償」というメリットは大きいのですが、Gaussian に比べると、利用者は少ないように思います
GAMESS に言及する場合には、３種類のなかで最も代表的なアメリカ版（GAMESS-US）について説明します
量子化学計算のコスト
量子化学計算プログラム
Gaussian
GAMESS
学術 / 商用：0 円
GUI アプリ
GaussView
学術：約 16 万円
商用：約 30 万円
計算用サーバー
約 17 万円 / 年
量子化学計算を比較的安価に実行できるスパコン環境
学術利用のみです
無償
産業利用枠もあります
成果を公開する場合、学術・産業利用のどちらも、１口（ = 1000 ノード時間）あたり 11万円
産業利用をメインにした計算環境です
実演：スパコンを用いたエタノールの量子化学計算
https://scrapbox.io/files/63fd3bf4fe79f6001cc36855.mov
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秒までの計算が可能
タイミングが悪いと混んでいる
この授業用
Username: 授業内でアナウンス
Password: 授業内でアナウンス
3分までの計算が可能
全体で8コアまで
分子を SMILES で指定するときには、メニューから...
Lookup → Import by → SMILES
実演：WebMO を用いたエチレンの量子化学計算
https://scrapbox.io/files/63f99217e5f976001c763dc8.mov
【演習課題３】WebMO を用いたエタノールの量子化学計算
WebMO （授業用）を用いて、エタノール（$ \mathsf{C_2H_5OH}）の量子化学計算（ GAMESS ） を実行してください。ただし、必ず、構造最適化（Geometry Optimization） を実行してください。計算手法には Hartree-Fock (RHF)、と基底関数には6-31G(d)を用いてください。計算が終了したら、得られた結果に基づいて、次に示す化学的性質を調べてみてください。 C-O 間の結合距離
Molecule Viewer
1.405 Å
全エネルギーの値
Overview → RHF Energy
-154.0757448250 Hartree
双極子モーメントの値
Overview → Dipole Moment
1.738 Debye
O 原子上の部分電荷の値
Partial Charges
-0.735
23-10-03 第３回
今回の内容
量子化学計算の基本的な手順
量子化学計算の代表的な計算方法
基本となるハートリー・フォック法
量子化学計算の基本的な手順
量子化学計算の基本的な手順
分子の状態（電荷とスピン多重度）を指定する
分子の初期配置を指定する
分子の計算方法と基底関数を指定する
分子の状態を指定する
電荷（Charge）
計算する系全体の総電荷を指定する
例えば...
１価の陽イオン（例：$ \mathsf{Na}^+）
電荷：１
１価の陰イオン（例：$ \mathsf{Cl}^-）
電荷：−１
一価の陽イオンと陰イオンが混在する場合（例：$ \mathsf{Na}^+ + \mathsf{Cl}^-）
電荷：０
スピン多重度（Multiplicity）
スピン多重度は、全スピン角運動を $ S とするとき、$ 2S+1で定義される
α電子とβ電子の電子数の差が $ 2Sになる
例えば...
不対電子を持たない化合物（例：基底状態の $ \mathsf{H_2}）
電子配置の例：↑↓
全スピン角運動量 $ S：０
スピン多重度：１（Singlet）
不対電子を２個もつ化合物（例：三重項状態の$ \mathsf{H_2}）
電子配置の例：↑↑
全スピン角運動量 $ S：１
スピン多重度：３（Triplet）
不対電子を１個もつ化合物（例：$ \mathsf{H_2^+}）
電子配置の例：↑
全スピン角運動量 $ S：1/2
スピン多重度：２（Doublet）
電荷とスピン多重度の例
table:電荷とスピン多重度
中性分子 カチオン アニオン ラジカル 励起三重項
電荷 0 1 -1 0 0
不対電子数 0 0 0 1 2
スピン多重度 Singlet Singlet Singlet Doublet Triplet
WebMOで分子の状態を指定する
https://i.gyazo.com/fe597d2dc6dd5a28b79fa97ddc5a509b.png
【演習課題１】分子の電荷とスピン多重度
次に示す分子の電荷とスピン多重度について、適切な値を答えてください。
基底状態のホルムアルデヒド $ \mathsf{CH_2O}
電荷：0
スピン多重度：Singlet
基底状態の酸素分子 $ \mathsf{O_2}
電荷：0
スピン多重度：Triplet
フェノールアニオン $ \mathsf{C_6H_5O^-}
電荷：-1
スピン多重度：Singlet
一酸化窒素 $ \mathsf{NO}
電荷：0
スピン多重度：Doublet
アンモニウムカチオン$ \mathsf{NH_4^+}
電荷：+1
スピン多重度：Singlet
分子の初期配置を指定する
量子化学計算では、必ず、分子の構造を入力する必要がある
入力する構造（初期構造）を出発点として、最も安定なエネルギーを持つ構造（最安定構造）を求める計算を「構造最適化」（Geometry Optimization）という
構造最適化の効率は、初期構造に大きく依存する
初期構造が適切ではないと、最安定構造（global minimum）ではなく、局所安定構造（local minimum）に収束してしまうことがある
例えば、ブタンの場合...
https://livedoor.blogimg.jp/commencement_wisdom/imgs/3/1/319eb5bb.jpg
分子の計算方法と基底関数を指定する
計算方法の種類 と 基底関数の種類 を スラッシュ（/） で挟んで表記することが多い
たとえば...
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）の組み合わせで、エネルギー計算を実行
WebMOで分子の計算方法と基底関数を指定する
https://i.gyazo.com/0fffdf97e906bc3d9f7412c20413606e.png
実演：WebMO を用いたブタンの構造最適化
https://i.gyazo.com/15ed58c53823a47f5b9a9eeaa6f4ae6f.mp4
量子化学計算の代表的な計算方法
代表的な計算方法
電子相関を考慮しない計算方法
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)
もちろん、ここにある$ \psi^0(\mathbf{r}_1)や$ \psi^0(\mathbf{r}_2)は、「これから求めたいもの」なので、その中身が分からない
➡ 拡張 Huckel 近似や半経験法など、精度の粗い方法をつかって、とりあえず求めておく
ステップ②：「電子２がつくる平均場」中にある電子１について、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)
ステップ①で「とりあえず」ということで決めた$ \psi^0(\mathbf{r}_2)をつかった計算なので、得られる答えも、精度的にはまだまだ不十分
ステップ③：「電子１がつくる平均場」中にある電子２について、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)
ステップ①で「とりあえず」ということで決めた$ \psi^0(\mathbf{r}_1)をつかった計算なので、得られる答えも、精度的にはまだまだ不十分
ステップ④：ステップ②と③で得られた軌道エネルギー $ \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
ステップ②の続きでは、アップデートした$ \psi^0(\mathbf{r}_2)をつかった計算なので、得られる答えは、１回目よりも改善される（エネルギーが下がる）
https://gyazo.com/f34b20c7a201c8c9aade6ffe65e3933b
上図は HeH+の場合
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 以下
化学反応を定量的に議論することが可能な精度
【演習課題２】Hartree-Fock 法
量子化学計算の基本となる「Hartree-Fock 法」の「メリット」と「デメリット」について、下記の QR コードのリンク先にある「投票」のページから投稿してください。
投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
【演習課題３】Hartree-Fock 法
Hartree-Fock (HF) 法に関する次の文章について、適切なものをすべて選んでください。
x HF 法は、「ある電子について、他の電子の作る平均化された場の（平均場）中で、他の電子とは独立に運動している」という描像で近似している。 _ HF 法は、電子同士の相互作用に由来する「電子相関」について、定量的な精度で十分に考慮することができる計算方法である。 x HF 法は、分子の立体構造や化学的性質について、多くの場合、十分な精度で予測できる。たとえば、全エネルギーの値についても、99%程度の精度で見積もることができる。 _ HF 法は、化学反応に伴う活性化エネルギーについて定量的な精度で議論したい場合でも、十分に信頼できる結果が得られる計算方法である。 _ HF 法は、複数の電子のふるまいについて、「ある電子は、他の電子の位置にあらわに依存しながら運動している」という多電子描像で記述している。 23-10-10 第４回
今回の内容
量子化学計算の代表的な計算方法
よく使われている密度汎関数法
よく使われている密度汎関数法
電子相関を考慮する計算方法
密度汎関数法
★ 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 を特徴付ける重要なもの
Kohn-Sham 法
KS 法を使わない DFT モデル
Thomas と Fermi は、ある系が局所的には一様密度の電子ガス（均一電子ガス）としてふるまうと近似した汎関数を提唱した
https://gyazo.com/fbc31958ac35ce0c52c1334a2b769b49
KS 法を使わない DFT モデル（T）の大きな欠点は...
運動エネルギー汎関数$ T[\rho] の記述が悪い
この問題に対して Kohn と Sham は...
電子密度$ \rho 使わず、1 電子波動関数（軌道）$ \phi を使って$ T[\rho] を表現するアイデアを提唱した
このアイデアを大まかに説明すると...
有効外部ポテンシャル $ v_s(\mathbf{r})のもとにある $ n個の相互作用がない仮想の電子系（Kohn-Sham 系）を考える
ただし$ v_s(\mathbf{r})は、相互作用を考慮したリアルな系で求められるものと等しい（仮想的な分子系）とする：
https://gyazo.com/9d050d7261c04608e767abdd31023288
Hartree-Fock 法とDFT法（Kohn-Sham法）の違い
KS 方程式とHF 方程式の違いは、３項目が「平均場$ V^\mathrm{HF}」か「外場有効ポテンシャル$ V^\mathrm{eff}」かのみ
つまり、KS 方程式は、HF 方程式と同程度の時間で計算できる
HF 方程式中の「平均場$ V^\mathrm{HF}」は「既知の関数」であるが、KS 方程式中の「交換相関汎関数$ V_\mathrm{xc}」は厳密な形が不明な「未知の関数」である
$ V_\mathrm{xc}の厳密な形が分からないので、人為的に仮定された汎関数を使う
この「人為的に仮定された」汎関数の質が、DFT の精度を決定することになる
【演習課題１】密度汎関数法
量子化学計算でよく用いられる「密度汎関数法」は、「Hartree-Fock 法」と比較したときに、どのような「メリット」と「デメリット」があると考えられるか。下記の QR コードのリンク先にある「投票」のページから投稿してください。
投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
密度汎関数法で用いる「汎関数」
汎関数とは？
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, ωB97X, ωB97X-D）
LC-BOP
分散力補正を含む汎関数
B3LYP-D3
ωB97X-D
ミネソタシリーズ（M06-2X, X11 など）
【演習課題２】水分子の計算
水分子（$ \mathsf{H_2O}）について、WebMO を使って、下記の表で指定される計算手法で構造最適化（Geometry Optimization）計算を実行してください。 table:h2o
学籍番号の末尾 Theory Functional Basis Set
0 RHF 6-31G(d)
1 RHF 3-21G
2 DFT BLYP 6-31G(d)
3 DFT BLYP 3-21G
4 DFT B3LYP 6-31G(d)
5 DFT B3LYP 3-21G
6 DFT PBE 6-31G(d)
7 DFT PBE 3-21G
8 DFT PBE0 6-31G(d)
9 DFT PBE0 3-21G
計算結果から、O-H 結合距離、および、H-O-H結合角度の値をそれぞれ読み取り、manaba の小テストにあるページに入力してください。
実験値：
O-H 結合距離：0.958 Å
H-O-H 結合角度：104.51 度
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 型に変えると、ほとんどの場合、計算精度が向上
【演習課題３】密度汎関数法
密度汎関数理論（DFT）法に関する次の文章について、適切なものをすべて選んでください。
_ DFT は、HF と同様に、波動関数を試行関数とする方法であり、その基礎となる Kohn-Sham 方程式は、シュレディンガー方程式を近似したものであると言える。 x DFT は、電子同士の相互作用に由来する「電子相関」について、Hartree-Fock 法と比較して、より定量的な精度で考慮することができる計算方法である。 x 一般的に、hybrid 型の GGA 汎関数は、非 hybrid 型の GGA 汎関数に比べて、様々な物性値を計算するときに、その精度が高いと期待できる。 _ DFT の汎関数として代表的なものは B3LYP であり、この汎関数は長距離補正を含むので、分子集合体系の分子間相互作用などについても、定量的な精度で正しく評価できる。 x ωB97X-D は、長距離補正だけではなく、分散力補正も含む hybrid 型の GGA 汎関数である。 23-10-17 第５回
今回の内容
量子化学計算の代表的な基底関数
量子化学計算の代表的な基底関数
基底関数とは？
任意の波動関数 $ \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://gyazo.com/ed83fb5d5c917d3e1f6f89baa7ad1716
基底関数には様々な種類があり、どのような分子を解析したいのか、何を調べたいと思っているのか、どのような計算方法を用いるのか、それぞれに応じて、適切に選択する必要がある
基底関数の個数が多いほど、電子分布を柔軟に表現できるので、計算精度が高くなる
計算精度が高くなればなるほど、計算コストも高くなってしまう
計算精度と計算コストのバランスを考えて、適切に選ぶ
【演習課題１】水分子の計算
水分子（$ \mathsf{H_2O}）について、WebMO を使って、下記の「割り当て」表で指定される計算手法でエネルギー（Molecular Energy）計算を実行してください。 計算結果から、「計算時間（Time）」を読み取り、下記の QR コード or このリンク先（削除）にあるページに、「計算方法、基底関数、計算時間」の３つを入力してください。
table:割り当て
学籍番号の末尾 Theory Functional Basis Set
0 RHF STO-3G
1 RHF 3-21G
2 RHF 6-31G(d)
3 RHF 6-311+G(d,p)
4 RHF cc-pVDZ
5 DFT B3LYP STO-3G
6 DFT B3LYP 3-21G
7 DFT B3LYP 6-31G(d)
8 DFT B3LYP 6-311+G(d,p)
9 DFT B3LYP cc-pVDZ
注意事項：
「構造最適化 (Geometry Optimization)」は、実行しないでください
Job Name のところに、必ず、自分の学籍番号を入力してください
https://gyazo.com/114bb9cb4c9163a71f418c84bd56a2a2
水素分子の場合
水素分子を表す波動関数 $ \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/5f46321dcacffafb09a3bacd2089e88e
ただし、数値計算の計算効率は悪い
yamnor.icon 図にあるように、原子の中心（$ r = 0）で値が不連続になるから、数値的に積分することが難しい
計算効率が悪いので、ふつうは、量子化学計算では使わない
ガウス型関数
量子化学計算では、スレーター型関数$ \exp(-\zeta r)の代わりに、$ \exp(-\alpha r^2)型の関数を基底関数として使う
$ \phi_a = \exp(-\alpha r^2)
このような関数を ガウス型関数 と呼ぶ
https://i.gyazo.com/2e928cc17ba5d8b6e8e1968aa8bdfc72.png
係数 $ \alphaは、関数の「広がりの程度」を表す変数
量子化学計算の前に、あらかじめ設定されている定数（計算の途中で変化させない）
ガウス型関数は、原子軌道に対する近似解なので、計算精度は高くはない
ただし、数値計算の計算効率が良い
yamnor.icon 上の図にあるように、原子の中心（$ r = 0）でも値が連続になるから、数値積分が簡単
計算精度の悪さは、拡がりの程度を変えたガウス型の基底関数をいくつか線形結合させることでカバーする
いくつかのガウス型関数を線形結合させたものを短縮型ガウス基底とよぶ
たとえば...
$ \phi_\mathsf{blue} = c_\mathsf{cyan} \, \chi_\mathsf{cyan} + c_\mathsf{pink}\, \chi_\mathsf{pink} + c_\mathsf{green}\, \chi_\mathsf{green}
$ = (1.33 e^{-0.59 r^2})_\mathsf{cyan} + (-0.57 e^{-1.28 r^2})_\mathsf{pink} + (0.04 e^{+0.06 r^2})_\mathsf{green}
https://i.gyazo.com/ec079f328422c5a601703e3c53390699.png
青い線で示す短縮型ガウス関数は、赤い線で示すスレーター型関数のかたちを、まあまあ良く再現している
基底関数の種類①
最小基底系（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) レベルの計算手法と組み合わせて用いる
【演習課題２】基底関数について
これまでの説明を聞いて、基底関数について「分かったこと」、「まだ分からないこと」は何か。
下記の QR コード or このリンク先（削除）にある「投票」のページから投稿してください。
投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
分極関数とは？
分子の複雑な構造を適切に記述したい場合には、各原子上に配置した基底関数に 分極関数（polarization function） を追加することで、電子分布が柔軟に変形する自由度を与えることができる
たとえば...
p 軌道をもつ水素以外の原子に対しては、d 型の関数を追加する
https://gyazo.com/40223e2a3c28a68693fe4ddf90a71bf5
球対称な s 軌道のみを持つ水素原子に対しては、p 型の関数を追加する
https://gyazo.com/a1cf4a2a749a46d7c53747177339d7ae
表記方法 / 指定方法
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
【演習課題３】基底関数について
基底関数に関する次の文章について、適切なものをすべて選んでください。
_ 量子化学計算では、ガウス型関数の代わりに、スレーター型関数を基底関数として使う。なぜならば、スレーター型関数は、計算精度が高いうえに、計算効率も良いからである。 x 最小基底系のひとつである STO-3G は、計算コストは低いが、計算精度は低いので、実験結果との定量的な比較などの実用的な解析に用いることは難しい。 x Pople 系基底関数のひとつである 6-311G は、一般に、同じ Pople 系基底関数である 3-21G よりも柔軟に電子のふるまいを記述できるため、計算精度が高い。 _ 高精度計算に適した基底関数のひとつである cc-pVDZ は、同じ関数系である cc-pVTZ よりも計算精度は高く、CCSD(T) レベルの計算手法と組み合わせて化学反応なども定量的に議論できる。 x Pople 系基底関数のひとつである 6-31G などは、固体表面や結晶などの周期性をもつ分子系の電子状態を解析しようとするときに使うことはあまりない。 x 水素分子（$ \mathrm{H_2}）についての量子化学計算をおこなう場合、基底関数として 6-31G と 6-31+G(d) のどちらを使っても、全エネルギーなどの計算結果は、まったく同じになる。 x フラーレン（$ \mathrm{C_{60}}）についての量子化学計算をおこなう場合、基底関数として 6-31+G(d) と 6-31++G(d,p) のどちらを使っても、全エネルギーなどの計算結果は、まったく同じになる。 _ 基底関数には様々な種類があるが、どのような場合でも、量子化学計算プログラムで用いることができる「最も高い精度の基底関数」を選べば良い。 _ 負電荷を帯びた分子は、中性分子と比較して電子分布が広がっているので、電子分布の広がりを記述するために、基底関数には 6-31G ではなく、分極関数を加えた 6-31G(d) を用いた方がよい。 23-10-24 第６回
今回の内容
量子化学計算に基づく構造最適化
量子化学計算に基づく構造最適化
分子の安定構造の探索
各原子核に働く力（エネルギー勾配 $ \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/bb347808513c0ddd4326458dee788cc6
安定構造を探索するシンプルな例
ある２原子分子について、原子間距離を $ r、最適構造での原子間距離を $ r_e、全エネルギー $ E(r)を次の式で近似できると考えてみる：
$ E(r) = \kappa (r - r_e)^2
$ \kappa：力の定数
原子間距離に対するエネルギー勾配の値は
$ \frac{dE(r)}{dr} = E'(r) = 2\kappa(r - r_e)
今、$ \kappa = 1、$ r_e = 1とすると、
$ E(r) = (r - 1)^2、$ E'(r) = 2(r - 1)
原子間距離$ r_\mathrm{old} = 2である初期構造を出発して、変位の大きさを決めるパラメータ$ \alpha = 0.1として、安定構造を探索するためのステップを進めると
$ r_\mathrm{new} = r_\mathrm{old} - \alpha E'(r_\mathrm{old}) = r_\mathrm{old} - 2 \alpha(r_\mathrm{old} - 1)
$ = 2 - 2 \times 0.1 \times 1 = 1.8
$ E(r_\mathrm{old}=2.0) = 1.00かつ$ E(r_\mathrm{new}=1.8) = 0.64なので、$ E_\mathrm{new} < E_\mathrm{old}となる
繰り返していくと...
https://scrapbox.io/files/6536c898421580001c18d5ac.png
https://scrapbox.io/files/6536c89cf47245001d2404cd.png
「カタチをちょっとずつ変化させる」を繰り返すことで、分子の立体構造を半自動的に最適化できる
実際には、多くの量子化学計算プログラムは、エネルギー勾配の方向に動かすという単純な方法（勾配降下法）ではなく、もっと効率的なアルゴリズムを使っている
Gaussian では、デフォルトでは、GEDIISというアルゴリズムを用いる
実演：Gaussian を用いた水分子の安定構造の探索
https://scrapbox.io/files/63f9a0bcd06e52001c9d4a71.mov
最適化される様子を確かめるため、HOH 結合角を 150° 程度にしたものを初期構造とした
【演習課題１】水分子の計算
水分子（$ \mathsf{H_2O}）について、WebMO を使って、下記の表で指定される計算手法で構造最適化（Geometry Optimization）計算を実行してください。 table:h2o
学籍番号の末尾 Theory Functional Basis Set
0 RHF 3-21G
1 RHF 6-31G(d)
2 DFT BLYP 3-21G
3 DFT BLYP 6-31G(d)
4 DFT B3LYP 3-21G
5 DFT B3LYP 6-31G(d)
6 DFT PBE 3-21G
7 DFT PBE 6-31G(d)
8, 9 DFT PBE0 3-21G
9 DFT PBE0 6-31G(d)
計算結果から、O-H 結合距離、および、H-O-H結合角度の値をそれぞれ読み取り、次に示す実験値との誤差を計算してください。
実験値：
O-H 結合距離：0.958 Å
H-O-H 結合角度：104.51 度
下記の QR コード or このリンク先にあるページに、「計算方法、基底関数、結合距離の誤差、結合角度の誤差」を入力してください。 注意事項：
Job Name のところに、必ず、自分の学籍番号を入力してください
https://gyazo.com/114bb9cb4c9163a71f418c84bd56a2a2
計算例：水分子の構造最適化
https://gyazo.com/c2e381a0aff6f0bfb93bac93b78c270f
分子の遷移状態の探索
ある化学反応について、反応速度を議論するときには、その反応に伴う遷移状態を調べる必要がある
たとえば、ある反応の速度定数 $ 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
【演習課題２】構造最適化について
これまでの説明を聞いて、構造最適化と遷移状態解析について「分かったこと」、「まだ分からないこと」は何か。
下記の QR コード or このリンク先にある「投票」のページから投稿してください。 投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
遷移状態で振動解析すると...
値が負（虚数）となる振動数がある
https://i.gyazo.com/049b4b6c5eb05f92272957fb07a7d691.png
安定構造では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
局所安定状態（下に凸の放物線）であり、エネルギーの座標に関する二次微分の値が常に正である
$ \frac{\partial^2 E}{\partial Q^2} > 0
https://gyazo.com/bb347808513c0ddd4326458dee788cc6
遷移状態では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
特定の変位方向に対して局所不安定状態（上に凸の放物線）であり、その変位方向にはエネルギーの座標に関する二次微分の値が負となる
$ \frac{\partial^2 E}{\partial Q^2} < 0
振動解析すると、１つの振動モードの振動数のみが負の値になる
【演習課題３】構造最適化
構造最適化に関する次の文章について、適切なものを全て選んでください。
_ 最適構造では、エネルギーの座標に関する一次微分（$ \frac{\partial E}{\partial Q}）は「常に負の値」をもつ。 _ 最適構造では、エネルギーの座標に関する二次微分（$ \frac{\partial^2 E}{\partial Q^2}）は「常にゼロ」になる。 x 遷移状態では、エネルギーの座標に関する一次微分（$ \frac{\partial E}{\partial Q}）は「常にゼロ」になる。 _ 遷移状態では、エネルギーの座標に関する二次微分（$ \frac{\partial^2 E}{\partial Q^2}）は「常に正の値」になる。 x 最適構造では、振動解析すると、すべての振動モードの振動数が正の値になる。 _ 遷移状態では、振動解析すると、すべての振動モードの振動数が負の値になる。 x 分子の安定構造を探索するとき、多くの量子化学計算プログラムでは、エネルギー勾配の方向に動かすという単純な方法ではなく、GEDIIS などの効率的なアルゴリズムを用いている。 x 分子の遷移状態を探索するとき、多くの量子化学計算プログラムでは、エネルギーの二次微分で得られる Hessian 行列に基づく方法を使うので、安定構造の探索よりも計算コストがかかる。 23-10-31 第７回
今回の内容
量子化学計算に基づく振動解析
量子化学計算に基づく振動解析
振動解析とは？
分子固有の振動数と赤外強度を理論的に計算する方法
実験的に測定された赤外スペクトルの解釈を補助する手段として活用されている
各振動モードに対するラマン強度を計算することも可能
振動状態に基づいて、分子の熱力学的な情報（自由エネルギーなど）を計算することが可能
➡ 振動状態を解析しないと、ある温度における自由エネルギーの大きさなど、分子の熱力学的な性質についての詳しい情報が得られない
最適化した構造が「安定点」か「不安定点（遷移状態）」かを見分けることにも使える
赤外吸収スペクトルとは？
水分子の赤外吸収スペクトル（実験結果）
https://gyazo.com/246c01937cec7d08b6c54b43d41e6f3e
分子振動とは？
水（$ \mathrm{H_2O}）の変角振動 / 対称伸縮振動 / 逆対称伸縮振動
https://gyazo.com/e4d9bb89f0f61f0533b9532d0e024287https://gyazo.com/145fab5f0795498d91045fd285f6b556https://gyazo.com/80f8eb416e3294f461be22f863b7cc58
分子振動の自由度
分子振動の自由度（振動モードの数）$ L_\mathrm{vib} は、$ N原子分子の場合...
直線分子：$ L_\mathrm{vib}=3N-5
非直線分子：$ L_\mathrm{vib}=3N-6
なぜならば...
$ L_\mathrm{vib}=L_\mathrm{all}-(L_\mathrm{trans}+L_\mathrm{rot})
分子の全自由度：$ L_\mathrm{all} = 3N
並進の自由度：$ L_\mathrm{rot} = 3
回転の自由度
直線分子：$ L_\mathrm{rot}=2
非直線分子：$ L_\mathrm{rot}=3
【演習課題１】振動自由度の数
水（H2O）、二酸化炭素（CO2）、メタノール（CH3OH）の振動自由度の数はいくつか？
水：3 × 3 - 6 = 3]
二酸化炭素：3 × 3 - 5 = 4]
メタノール：3 × 6 - 6 = 12]
分子振動を表す式（古典力学の場合）
原子１と２から成る２原子分子を考える
https://gyazo.com/732ce74906d08bff901aaa86ceaa1f19
原子間の化学結合をバネとみなして、力の定数を$ fとして、フックの法則が成り立つと仮定する
この場合のニュートン方程式は...
$ \mu\frac{d^2\Delta r}{dt^2}=-f\Delta r
$ \muは換算質量と呼ばれ、定義は...
$ \mu=\frac{m_1 m_2}{m_1 + m_2}
上式は単振動を表す運動方程式であり、その解は...
$ \Delta r = A \cos(\omega t)
$ \omega = \sqrt{\frac{f}{\mu}}
$ \omegaは角振動数と呼ばれ、振動数$ \nuとの関係は...
$ \omega = 2\pi v
したがって、振動数$ \nuで表すと...
$ \nu=\frac{1}{2\pi}\sqrt{\frac{f}{\mu}}
分子振動を表す式（量子力学の場合）
分子振動を調和振動子であるとして近似（調和近似）して扱う 振動状態の波動関数を$ \Psi_\mathrm{vib}とすると、対応するシュレディンガー方程式は...
$ \hat{H}_\mathrm{vib}\psi_\mathrm{vib}=\left( -\frac{\hbar^2}{2\mu}\frac{d^2}{d(\Delta r)^2} +\frac{1}{2}\mu\omega^2(\Delta r)^2\right) \psi_\mathrm{vib}=E_\mathrm{vib}\psi_\mathrm{vib}
シュレディンガー方程式を解くと、エネルギーは...
$ E_\mathrm{vib}=hv\left(n+\frac{1}{2}\right)
ここで $ n = 0, 1, 2, \cdotsは、振動量子数と呼ばれる
振動量子数$ nに対する、振動波動関数 $ \Psi_nとエネルギー$ E_nは...
https://gyazo.com/86b99039bd5170e0d4a80d90c2a99097
振動解析の計算手順
対象とする分子について、構造最適化を実行（安定構造を探索）する
構造最適化を実行したときと同じ量子化学計算レベルで、振動解析を実行する
実演：Gaussian を用いた水分子の振動解析
https://scrapbox.io/files/63f9de9406a787001b9cc00a.mov
【演習課題１】水分子の計算
水分子（$ \mathsf{H_2O}）について、WebMO を使って、下記の表で指定される計算手法で振動解析（Vibrational Frequencies）計算を実行してください。その際、必ず、振動解析の前に「構造最適化（Geometry Optimization）」をおこなってください。 table:h2o
学籍番号の末尾 Theory Functional Basis Set
0 RHF 6-31G(d)
1 RHF 3-21G
2 DFT BLYP 6-31G(d)
3 DFT BLYP 3-21G
4 DFT B3LYP 6-31G(d)
5 DFT B3LYP 3-21G
6 DFT PBE 6-31G(d)
7 DFT PBE 3-21G
8 DFT PBE0 6-31G(d)
9 DFT PBE0 3-21G
計算結果から、変角振動、対称伸縮振動、逆対称伸縮振動の値をそれぞれ読み取り、次に示す実験値との誤差を計算してください。
実験値：
変角振動：1595 cm-1
対称伸縮振動：3657 cm-1
逆対称伸縮振動：3756 cm-1
下記の QR コード or このリンク先にあるページに、「計算方法、基底関数、変角振動の誤差、対称伸縮振動の誤差、逆対称伸縮振動の誤差」を入力してください。
リンク削除
注意事項：
Job Name のところに、必ず、自分の学籍番号を入力してください
https://gyazo.com/114bb9cb4c9163a71f418c84bd56a2a2
計算例：水分子の振動解析
https://gyazo.com/288677666289ced8a6fec15061e4bcb8
振動解析に基づいた熱力学量の見積もり
分子の振動状態に基づいて、電子状態のエネルギーに対して、 熱力学量を「近似的に」算出するための補正をおこなうことができる
振動数が低い振動モードほど、熱力学補正に対する寄与が大きい
「近似」なので、大きく構造が揺らぐ柔軟な分子については、誤差が大きくなる場合がある
求めることができるのは...
ゼロ点エネルギー（Zero Point Energy, ZPE）振動補正
$ E_0 = E + \mathrm{ZPE}
熱力学補正
$ E_\mathsf{total} = E_0 + E_\mathsf{並進} + E_\mathsf{回転} + E_\mathsf{振動}
エンタルピー補正
$ H = E_\mathsf{total} + k_\mathrm{B} T
Gibbs 自由エネルギー補正
$ G = H - T S_\mathsf{total}
$ S_\mathsf{total} = S_\mathsf{並進} + S_\mathsf{回転} + S_\mathsf{振動} + S_\mathsf{電子}
熱力学補正の具体例
$ \mathsf{C_2H_5 + H_2 \rightarrow C_2H_6 + H}
https://gyazo.com/eb01b32e0248d4549310ae074b5f5068
エネルギー変化（熱力学補正をしない場合）
$ \Delta E = E_\mathsf{生成物} - E_\mathsf{反応物}
4.86 kcal/mol
反応エンタルピー
$ \Delta_r H^\circ(298\mathsf{K}) = H_\mathsf{生成物} - H_\mathsf{反応物}
8.08 kcal/mol
Gibbs 自由エネルギー変化
$ \Delta_r G^\circ(298\mathsf{K}) = G_\mathsf{生成物} - G_\mathsf{反応物}
11.18 kcal/mol
分子振動の非調和性
高い精度（MP2/cc-pVTZ）で計算した場合
https://gyazo.com/e1d212dc0d548ac021fe5ffb7701964a
MP2摂動法を用いることで、実験値にかなり近くなる
高波数の振動モード（対称伸縮＆逆対称伸縮）では、誤差が大きい
分子振動を調和振動子だと見なす近似（調和近似）による誤差（非調和性）の影響だと考えられる https://i.gyazo.com/2d4bfe9a683c74a99f3f3098a5bebc3f.png
非調和性の問題を手軽に解決する場合には、振動数に適当な数字（スケーリング因子）を掛けて補正する
事例：
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
【応用事例】結晶多形を振動スペクトルで見分ける
HPBI という分子の結晶多形（左：α型 / 右：β型）
分子内の水素結合のある・なしで、3000〜3500 cm-1 のバンド形状が大きく異なることが分かりました
https://gyazo.com/707a1f02e3588843dd5b74e8d3408537 https://gyazo.com/4fd67a5a3962e915c0cd10e1dcdddbbd
https://gyazo.com/7c73f8dab6b4f621a2a8ff043a2cd130https://gyazo.com/5d4566bc2ab7d177e39c89c37edc099d
https://gyazo.com/2856c129d532d28ecf76bc2969f3f604 https://gyazo.com/4f407938cb8724110cf4dfd9ae03451a
山本@山本研, 論文準備中
【演習課題３】振動解析について
これまでの説明を聞いて、振動解析について「分かったこと」、「まだ分からないこと」は何か。
下記の QR コード or このリンク先にある「投票」のページから投稿してください。
リンク削除
投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
23-11-07 第８回
今回の内容
量子化学計算に基づく化学反応の解析
量子化学計算に基づく化学反応の解析
量子化学計算によって実験結果の予測、実験の手助けとなる計算を行うことは？
化学反応にともなう活性化エネルギー（反応障壁の大きさ）を見積もることで、「反応の速度」（反応速度定数）を理論的に見積もることができます
反応速度定数の見積もりには、「アレニウス則」や「遷移状態理論」などを用います
また、反応物と生成物のエネルギー差（反応エンタルピー）から、ある温度における「反応物と生成物の分布比」を理論的に見積もることができます
分布比は「ボルツマン分布」に従います
分子の遷移状態
ある化学反応について、反応速度を議論するときには、その反応に伴う遷移状態を調べる必要がある
たとえば、ある反応の速度定数 $ kについて、始状態と遷移状態とのエネルギー差、つまり、活性化エネルギー $ \Delta G^{\ddagger}が分かっていれば、Eyring-Polanyi 式から...
$ k = \chi \frac{k_\mathrm{B}T}{h}\exp\left(-\frac{\Delta G^{\ddagger}}{RT}\right)
$ \chi：遷移状態において、反応座標に沿った分子振動が一往復したときに反応する確率（透過係数）
遷移状態に基づく解析の具体例
$ \mathsf{FH + Cl \rightarrow F + H Cl}
水素（$ \mathsf{H}）の場合
$ \Delta G^{\ddagger}= 17.26 kcal/mol
$ k(298\mathsf{K}) = \frac{k_\mathrm{B}T}{h}\exp\left(-\frac{\Delta G^{\ddagger}}{RT}\right) = 1.38 $ \mathsf{s}^{-1}
透過係数$ \chiは 1 として計算
重水素（$ \mathsf{D}）に置換した場合
$ \Delta G^{\ddagger}= 18.86 kcal/mol
$ k(298\mathsf{K}) = 0.0928 $ \mathsf{s}^{-1}
➡ 重くなることで、反応速度が遅くなる
【実演】Gaussian を用いた Diels-Alder 反応の解析
https://gyazo.com/6245b572e337bc0273b0b260af89a429
https://i.gyazo.com/a53e519a0a2b7c54971c9b3bd1ef1871.png
【応用事例】食品成分（リコピン 🍅）の trans → cis 異性化を促進したい
https://gyazo.com/84c251e60a9a9d9eeb71ddc2a807c97d
https://gyazo.com/6b1cd074ab2225a2ac3d0b11b7486e23
宮川@山本研, 分子科学討論会 (2023)
【応用事例】色素材料の円偏向性を制御したい
凝集誘起発光特性をもつ三脚巴状分子の「らせん反転運動」の制御
https://gyazo.com/e23ca93fd21e4bb92b9e7671dfa9d84c
https://gyazo.com/599f6514b065d5694402512301bd2e68
大場@山本研, 分子科学討論会 (2023)
【応用事例】酵素反応のメカニズムを解析したい
https://gyazo.com/d1e23cee45d8ee32d7f2e9580e97fcbb https://gyazo.com/c24ebd0b6dd52bbfe79bfffd6d3d2cd4
山本＠山本研, Life (2022)
遷移状態の探索
多くの量子化学計算プログラムには、遷移状態を探索する機能が備わっている
単純に「構造最適化の逆」（$ x_\mathrm{new} = x_\mathrm{old} + \alpha \frac{\partial E}{\partial Q}）をやっても、うまくいかない
https://gyazo.com/bb347808513c0ddd4326458dee788cc6
エネルギーの二次微分（$ \frac{\partial^2 E}{\partial Q^2}）を要素にもつ Hessian 行列 の固有値（力の定数 or 曲率：$ \chi）を求めて...
Hessian 行列の負の固有値を持つ固有ベクトルの方向に少しずつ動かすという方法がある
→ Eigenvector Following 法
調和近似（$ E = \kappa Q^2）からの歪みが大きな方向に少しずつ動かすという方法がある
Hessian 行列を計算するためのコストがまあまあ大きい
➡ 構造最適化よりも計算コスト（時間）がかかる
反応の場などの前提条件は、事前に設定したうえで遷移状態を探索するのでしょうか？
あらかじめ、「遷移状態に近いのはどのような構造か？」について、大まかな知識や前提が必要となります
実演：Gaussian を用いたブタジエンの cis/trans 異性化反応の解析
https://scrapbox.io/files/63faad6c14c7b1001bf5afac.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
【演習課題Ｘ】ブタジエンの cis/trans 異性化反応の解析
ブタジエンの cis/trans 異性化反応について、量子化学計算で得られた結果に基づいて、以下の問いに答えてみてください。量子化学計算を実行する際には、基底関数は 6-31G(d) 、計算方法は密度汎関数理論（DFT）とし、汎関数はB3LYPとします。
問１
ブタジエンの trans 体のエネルギー $ \epsilon_\mathrm{trans}を Hartree 単位で答えてください
$ \epsilon_\mathrm{trans}=-155.992139 Hartree
問２
ブタジエンの cis 体のエネルギー $ \epsilon_\mathrm{cis}を Hartree 単位で答えてください
$ \epsilon_\mathrm{cis}=-155.986486 Hartree
問３
ブタジエンの cis/trans 異性化反応について、遷移状態におけるエネルギー $ \epsilon_\mathrm{TS}を Hartree 単位で答えてください
$ \epsilon_\mathrm{TS}=-155.980090 Hartree
【演習課題１】ブタジエンの cis/trans 異性化反応の解析
ブタジエンの cis/trans 異性化反応について、量子化学計算で得られた結果に基づいて、以下の問いに答えてみてください。ただし、気体定数の値は$ R = 1.987\,206 \; \mathrm{cal} \; \mathrm{K}^{-1} \; \mathrm{mol}^{-1}、温度は$ T = 300 \; \mathrm{K}、頻度因子の値は$ A = 1.0 \times 10^{14} \; \mathrm{sec}^{-1}とします。また、エネルギーの単位である Hartree については、1 Hartree $ = 627.5095 \; \mathrm{kcal} \; \mathrm{mol}^{-1}となります。
問１
ブタジエンの trans 体と cis 体の存在比$ \frac{P_\mathrm{trans}}{P_\mathrm{cis}}の値を答えてください
$ \Delta E_\mathrm{trans-cis}= -3.547 $ \mathrm{kcal \; mol^{-1}}
$ \frac{P_\mathrm{trans}}{P_\mathrm{cis}} = \exp\left(-\frac{\Delta E_\mathrm{trans-cis}}{RT}\right)
$ = \exp\left(-\frac{\Delta E_\mathrm{trans-cis} \times 10^{3}\;\mathrm{cal\;mol^{-1}}}{1.987\,206 \; \mathrm{cal\; K^{-1} \; mol^{-1}}\times 300 \; \mathrm{K}}\right)
$ =384
問２
ブタジエンの trans → cis 異性化反応に対する反応速度定数$ k_{\mathrm{trans}\rightarrow\mathrm{cis}}を$ \mathrm{sec}^{-1}単位で答えてください
$ \Delta E_\mathrm{TS-trans} = 7.561 $ \mathrm{kcal \; mol^{-1}}
$ k_{\mathrm{trans}\rightarrow\mathrm{cis}} = A \exp\left(-\frac{\Delta E_\mathrm{TS-trans}}{RT}\right)
$ = 1.0 \times 10^{14} \; \mathrm{sec}^{-1} \; \times \exp\left(-\frac{\Delta E_\mathrm{TS-trans} \times 10^{3}\;\mathrm{cal\;mol^{-1}}}{1.987\,206 \; \mathrm{cal\; K^{-1} \; mol^{-1}}\times 300 \; \mathrm{K}}\right)
$ = 3.105 $ \times 10^{8} $ \mathrm{sec^{-1}}
問３
ブタジエンの cis → trans 異性化反応に対する反応速度定数$ k_{\mathrm{cis}\rightarrow\mathrm{trans}}を$ \mathrm{sec}^{-1}単位で答えてください
$ \Delta E_\mathrm{TS-cis} = 4.014 $ \mathrm{kcal \; mol^{-1}}
$ k_{\mathrm{trans}\rightarrow\mathrm{cis}} = A \exp\left(-\frac{\Delta E_\mathrm{TS-cis}}{RT}\right)
$ = 1.0 \times 10^{14} \; \mathrm{sec}^{-1} \; \times \exp\left(-\frac{\Delta E_\mathrm{TS-cis} \times 10^{3}\;\mathrm{cal\;mol^{-1}}}{1.987\,206 \; \mathrm{cal\; K^{-1} \; mol^{-1}}\times 300 \; \mathrm{K}}\right)
$ = 1.192$ \times 10^{11} $ \mathrm{sec^{-1}}
遷移状態で振動解析すると...
値が負（虚数）となる振動数がある
https://i.gyazo.com/049b4b6c5eb05f92272957fb07a7d691.png
安定構造では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
局所安定状態（下に凸の放物線）であり、エネルギーの座標に関する二次微分の値が常に正である
$ \frac{\partial^2 E}{\partial Q^2} > 0
https://gyazo.com/bb347808513c0ddd4326458dee788cc6
遷移状態では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
特定の変位方向に対して局所不安定状態（上に凸の放物線）であり、その変位方向にはエネルギーの座標に関する二次微分の値が負となる
$ \frac{\partial^2 E}{\partial Q^2} < 0
振動解析すると、１つの振動モードの振動数のみが負の値になる
遷移状態を見つけるのが難しい場合には？
Nudged Elastic Band (NEB) 法を用いて、最小エネルギー経路を探索する
https://gyazo.com/e43f8e5878a51ef4e97d8f9106c24fc4https://gyazo.com/d42120239cf54afc31cf6068fccabd06
Gaussian で使える？
現在のバージョンには実装されていない
HPC システムズ（株）が開発・販売している Reaction Plus（学術：80 万円、商用：200 万円）というソフトウェアと Gaussian を組み合わせると NEB 計算ができる 【演習課題２】化学反応の解析について
これまでの説明を聞いて、量子化学計算に基づく化学反応の解析について「分かったこと」、「まだ分からないこと」は何か。
下記の QR コード or このリンク先にある「投票」のページから投稿してください。
リンクは削除
投稿するときには、複数の答えがある場合にもまとめて、１回で投稿してください。
この演習は、成績評価の対象となります。投票者の「名前」を、匿名から「学籍番号」に変更してください。
23-11-14 第９回
今回の内容
量子化学計算に基づく電子励起状態の解析
量子化学計算に基づく溶媒効果の解析
量子化学計算に基づく電子励起状態の解析
🟧【事前質問】化合物の反応性について、熱や光など、化合物の特性以外の要因の影響は見ることは？
「光」については、化合物を光で励起した後の状態（励起状態）について、量子化学計算を用いて、その影響を調べることが可能です
たとえば...
化合物が吸収・発光する光の波長 / 強度
光励起したあとの化合物の構造変化・反応過程
光の吸収とは？
光のエネルギーが分子を基底状態から励起状態に変えるエネルギーとして消費されるプロセス
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}）状態
つまり、量子化学計算に取り組むときに、分子の立体構造は、励起状態で最適化する必要がある
【応用事例】高効率な有機ＥＬを設計したい
分子軌道のカタチを調べることでロジカルな設計が可能となった
➡ HOMO と LUMO の重なりを少なくすることで、$ \mathsf{T_1}→ $ \mathsf{S_1}の項間交差（ISC）を促進する
https://gyazo.com/bb987be1f6b48f506a80c7d00a616cf0
安達＠九大, et al., Nature Comm., Vol. 6, p. 8476 (2015)
【応用事例】環境応答性をもつ色素材料を設計したい
https://gyazo.com/ab7a6e54a31a4748fa357b47215a0d7a
https://gyazo.com/647ee04cd1561997bf49c4fcc4122f73
古川@山本研, 分子科学討論会 (2023)
【応用事例】古代ホタルは何色に光るのか？
https://gyazo.com/df5049e3be79d42cc3f417d4441e5357
https://gyazo.com/33fec10dfe56cbccc1bfd9b1ad889d94
尾保手@山本研, 分子科学討論会 (2023)
励起状態の計算方法
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とは？
時間依存 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
【演習課題１】ホルムアルデヒドの励起状態の計算
ホルムアルデヒド（$ \mathsf{H_2CO}）について、WebMO を使って、下記の表で指定される計算手法で励起エネルギー Energy (Excited State)計算を実行してください。その際、必ず、振動解析の前に「構造最適化（Geometry Optimization）」をおこなってください。 table:h2o
学籍番号の末尾 Theory Functional Basis Set
0, 5 RHF 6-31G(d)
1, 6 DFT BLYP 6-31G(d)
2, 7 DFT B3LYP 6-31G(d)
3, 8 DFT PBE 6-31G(d)
4, 9 DFT PBE0 6-31G(d)
この分子の吸収スペクトルで観測される電子遷移のうち、最も光吸収の程度が大きいと予測される遷移に対応する「吸収波長の値（nm 単位）」と「振動子強度の値」を教えてください。
下記の QR コードにあるページに、下記の例を参考にして、「計算方法（例：RHF）、吸収波長（例：123 nm）、振動子強度（例：0.456）」を入力してください。
code:参考
RHF ： 123 nm （0.456）
リンクは削除
注意事項：
Job Name のところに、必ず、自分の学籍番号を入力してください
https://gyazo.com/114bb9cb4c9163a71f418c84bd56a2a2
量子化学計算に基づく溶媒効果の解析
🟧【事前質問】
反応溶媒を選択する際に何か指標となるような計算ができないか？
量子化学計算で目的物の選択性を向上させるためのヒントは？
溶液中での溶解度・溶解速度や、構造の取り方（イオン化の有無 / 単一分子 o r複合体で溶解しているかなど）の予測は？
溶媒和とは
溶質分子が、静電気力や水素結合などによって結びつき取り囲むことで溶質が溶媒中に拡散する現象
https://scrapbox.io/files/653ad5ae76d858001c30094a.png
溶媒効果とは？
化学において、反応性もしくは分子の会合に対して溶媒が及ぼす影響を指す
溶媒は溶解度、安定性、反応速度に影響を及ぼすため、適切な溶媒を選択することにより、化学反応を熱力学的・速度論的に制御できる
溶媒効果の具体例
アセチルアセトンの互変異性に対する平衡定数 $ K_\mathrm{keto-enol} = \mathrm{\frac{[enol]}{[keto]}}
https://scrapbox.io/files/653ad6108e5379001c66ebe4.png
table:acetylacetone
溶媒の種類 keto-enol 平衡定数
水 0.23
ジクロロメタン 4.2
エタノール 5.8
テトラヒドロフラン 7.2
気相 11.7
ベンゼン 14.7
シクロヘキサン 42
平衡定数とギブス自由エネルギー変化
ある可逆的な反応過程 $ A について、対応する平衡定数 $ Kとギブス自由エネルギー変化 $ \Delta Gの関係
$ \Delta G = -RT \ln K
$ R：気体定数（$ = k_\mathrm{B} N_\mathrm{A} = 8.314\,462 \; \mathrm{J \; K^{-1} \; mol^{-1}}）
$ T：温度（300 K）
【演習課題２】アセチルアセトンの互変異反応に対する溶媒効果
アセチルアセトンの互変異反応 $ \mathrm{keto \rightarrow enol} について、ある溶媒中におけるギブス自由エネルギー変化を測定したところ、 $ \Delta G = -6.70\; \mathrm{kcal\;mol^{-1}} でした。この場合、keto 体と enol 体の存在比はどのようになると予測できますか。
「manaba の小テスト」から回答してください。
量子化学計算で溶媒効果を考慮するには
溶媒分子を露わに考慮する（explicit solvation model）
https://scrapbox.io/files/653ad724f69f4f001b782be2.png
溶媒分子を分極する連続誘電体として扱う（implicit solvation model）
https://scrapbox.io/files/653ad7607fe76c001b737387.png
左：Debye-Onsager モデル
右：連続誘電体モデル
【演習課題３】メタノールの溶媒和エネルギー
メタノール（$ \mathsf{CH3OH}）について、WebMO を使って、下記の表で指定される種類の溶媒で溶媒和エネルギー（$ \Delta E_\mathrm{sol} = E_\mathrm{sol} - E_\mathrm{gas}）を計算してください。その際、必ず、エネルギー計算の前に「構造最適化（Geometry Optimization）」をおこなってください。計算方法は、密度汎関数理論（DFT）法、汎関数にはB3LYPを用いてください。基底関数は 6-31G(d) とします。また、1 Hartree = 627.5095 kcal/mol です。 table:CH3OH
学籍番号の末尾 溶媒の種類
0, 1, 2 Water
3, 4, 5 Acetoritrile
6, 7, 8, 9 Cyclohexane
下記の QR コードにあるページに、下記の例を参考にして、「溶媒の種類（例：Water）、溶媒和エネルギー（例：12.3 kcal/mol）」を入力してください。
code:参考
Water ： 12.3 kcal/mol
リンクは削除
アンケート
この授業では、授業の資料を Scrapbox でシェアしたり、WebMO を使って量子化学計算に取り組んでみたり、Slido をつかって演習課題の結果をシェアしてみたりと、「教員がしゃべるだけ」の授業からアップデートするための工夫をいろいろとやってみました。
このような授業の進め方は、みなさんにとって、役に立ったのでしょうか。
下記のアンケートから、みなさんの意見を教えてください。
みなさんからの意見を参考にして、より良い授業が実践できるように、いろいろと試行錯誤したいと思っています。
リンクは削除
https://gyazo.com/92d999fe1c25b85098eea4a2cb38a5e0