240226
https://gyazo.com/70e421419d3c18011ba6506d5ec80f8f
資料へのリンク
https://gyazo.com/5eee735b71f173d7fdd917f5fb05c93d
概要
日程： 2024 年 02 月 26 日（月）
会場： オンライン（Teams）
主催： 日本テクノセンター
目次
はじめに
量子化学計算の概要
量子化学計算の代表的なソフトウェア
量子化学計算の代表的な計算方法
量子化学計算の代表的な基底関数
量子化学計算に基づく構造最適化
量子化学計算に基づく振動解析
量子化学計算に基づく化学反応の解析
量子化学計算に基づく電子励起状態の解析
量子化学計算に基づく溶媒効果の解析
質疑応答
さいごに
講義① 10:30 〜 12:00
はじめに
学びの手引き
この資料について
資料は、講座が始まるまでに 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 の特徴と基本的な計算手順について、説明できる
量子化学の知識や技術に基づいて、おもに有機化合物の研究開発に応用するための指針について、説明できる
量子化学計算を活用して実験前の事前検証をおこない、研究・開発に必要となる時間的・資源的コストを効率化するための方針について、説明できる
この講座のタイムスケジュール
table:schedule
開始 終了 時間 内容
10:30 12:00 90分 講義①
12:00 13:00 60分 昼食
13:00 14:50 110分 講義②
14:50 15:10 20分 休憩
15:10 17:00 110分 講義③
17:00 17:30 30分 質疑
17: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 は無償配布
🟧【Ｑ】化合物が固体状態でも反応予測は可能なのでしょうか？
化合物が固体状態の場合、もし量子化学計算で扱おうとする場合には、一部を切り出したモデル系として扱うことが多いです（「クラスターモデル」と呼ぶことがあります）
固体であることの物性（バンド構造など）が重要となる場合には、「固体物理」の分野でよく用いられている、Quantum Espresso などの電子状態計算プログラムをつかったほうが、解析が進むかも知れません また、動的なふるまいに興味があったり、構造の揺らぎなどを詳しく知りたいときには、LAMMPS などの分子動力学シミュレーションなどを利用する必要があるかも知れません 🟧【Ｑ】金属吸着やPET等の各材質への吸着性を数値化して算出することは可能でしょうか？
金属表面への吸着を調べる場合には、上記でも少しご紹介した、クラスターモデルというアプローチで、数十個程度の金属でモデル化した表面に分子を置いて、その相互作用を調べることができます
吸着反応では、熱的な拡散過程など、動的なふるまいも重要になってきますので、量子化学計算よりも、分子動力学シミュレーションを用いることが多いかと思います
分子動力学シミュレーションについて、無料で使えるものでしたら、固体表面だと LAMMPS、多孔質材料だと RASPAなどが良いかと思います 【応用事例】金属有機構造体 MOF でゲスト分子がどのように拡散するのか
LAMMPS を用いた分子動力学シミュレーション
https://gyazo.com/c5d457492844fea35ecca6b1152e74c5
井上＠山本研, 日本コンピュータ化学会（2023）
【応用事例】インフルエンザウイルスの薬剤耐性化はどのようにおこるのか？
Amber を用いた分子動力学シミュレーション
https://gyazo.com/28406ac07980e525221bf2a2e861c2fe
Mohini＠山本研 et al, PeerJ (2021)
力試しテスト
次に示すような「新しい薬を開発するケース」について、量子化学計算を用いることが適切だと思われるものをすべて選んでください。
x 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子を合成するための準備として、合成経路の素反応のひとつについて、対応する遷移状態のエネルギーを知りたい。 x 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子について、分子内でどのような電荷分布になっているのか知りたい。 _ 新しい薬を設計するときに、１００万個の化合物データベースの中から、薬物候補化合物としてふさわしい分子を１０００個選びたい。 x 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子の赤外吸収スペクトルについて、それぞれのピークがどのような振動モードに由来するのかを調べたい。 _ 新しい薬を設計するときに、薬物候補化合物として有力だと思われる幾つかの分子について、どのような結晶構造になっているのかを予測したい。 量子化学計算の代表的なソフトウェア
業界標準の 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:ethylene.gjf
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 万円
学術 / 商用：0 円
学術：1,200 USD（約 16 万円）
商用：3,000 USD（約 41 万円）
計算用サーバー（クラウドサービス）
約 17 万円 / 年
約 1.2 千円 / 日
コストを掛けずに量子化学計算をはじめる
ちょっと試す
GAMESS + WebMO Basic + AWS
約 1.2 千円 / 日
もう少し本格的に
GAMESS + WebMO Pro (商用) + Sakura VPS
約 17 万円 / 年 + 約 41 万円
昼休憩 12:00 〜 13:00
ご質問①
GAMESSとGaussianの差について、性能面では異なる点はなにがありますでしょうか？
両者について、たとえば、分子の安定な構造を探索する「構造最適化」や分子の振動状態を解析する「振動解析」など、量子化学計算の基本的な手法については、どちらでも実行することができます
しかし、同じ計算手法で計算したのに、少し異なる値になってしまうこともあります。これは、繰り返しが必要な計算などをどの程度の閾値で打ち切るかなどの設定が異なるためです。
性能面では、「計算を実行するときのスピード」や「扱いやすさ」については、コストを掛けて開発されている Gaussian の方に軍配が上がるような気がします。
同じ「構造最適化」でも、Gaussian では「収束性の良いアルゴリズム」が実装されていることもあり、GAMESSよりも「計算が早く終わる」ことが多いように思います
また、ユーザー数としては Gaussian の方が多いように思いますので、たとえば、「論文に載っている方法と同じ手法をつかいたい」というような場合は、Gaussian の方が便利かも知れません
一方で、GAMESS には、Gaussian にないような「新しい計算手法」が積極的に導入される傾向にあります
たとえば、タンパク質などの大規模系を小さな断片（フラグメント）に分割して計算するフラグメント分子軌道法は、GAMESS でしか利用できません
ご質問②
結晶構造については、他手法で計算してからGaussianで計算が良いとのことですが、可能であれば他手法を具体的にご教示いただけないでしょうか？
Gaussianは、結晶構造など、周期性をもつ系も扱うことはできますが、「得意」ではありません
効率的に計算する場合には、他のプログラムを使う方が良いかもしれません
Gaussian → 原子基底
周期系を扱うときは → 平面波基底
分子性結晶の場合
分子性結晶の「構造探索」については、例えば、分子力場を用いて効率的に探索する Conflex というプログラム（有料）があります 分子性結晶の「電子状態解析」などについては、例えば、量子化学計算を用いて精度良く計算する Crystalというプログラム（有料）があります 金属結晶の場合
金属結晶の構造変化などは、分子力場を用いる LAMMPS （無料）が使えそうです 以上のようなものに加えて、最近では、第一原理計算のデータを深層学習したニューラルネットワークポテンシャル（NNP）力場を使ったプログラム（Matlantis や NeuralMD など）も登場してきました 第一原理計算に基づくものよりも５〜６倍のスピードで、第一原理計算レベルの精度で解析ができる場合があります
講義② 13:00 〜 14:50
量子化学計算の代表的な計算方法
代表的な計算方法
電子相関を考慮しない計算方法
Hartree-Fock（HF）法
電子相関を考慮する計算方法
密度汎関数法
★ Møller-Plesset 摂動法（例：MP2, MP4）
★ 結合クラスター法（例：CCSD, CCSD(T)）
★ 配置間相互作用（例：CISD, Full CI）法
実演：Gaussian で利用できる様々な計算方法
https://gyazo.com/717e9133b8a207143ca1e1288dd1bb3c
代表的な計算方法の精度とコスト
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(\mathbf{r}_1)とエネルギー $ \epsilon_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(\mathbf{r}_2)とエネルギー $ \epsilon_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 (HF) 法に関する次の文章について、適切なものをすべて選んでください。
x HF 法は、「ある電子について、他の電子の作る平均化された場の（平均場）中で、他の電子とは独立に運動している」という描像で近似している。 _ HF 法は、電子同士の相互作用に由来する「電子相関」について、定量的な精度で十分に考慮することができる計算方法である。 x HF 法は、分子の立体構造や化学的性質について、多くの場合、十分な精度で予測できる。たとえば、全エネルギーの値についても、99%程度の精度で見積もることができる。 _ 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 を特徴付ける重要なもの
このアイデアを大まかに説明すると...
外部有効ポテンシャル $ V^\mathrm{eff}のもとにある $ n個の相互作用がない仮想の電子系（Kohn-Sham 系）を考える
ただし$ V^\mathrm{eff}は、相互作用を考慮したリアルな系で求められるものと等しい（仮想的な分子系）とする：
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 の精度を決定することになる
密度汎関数法で用いる「汎関数」
汎関数とは？
DFT の計算精度と信頼性は、電子密度と系のエネルギーを関係付ける汎関数の種類に依存する
適切な汎関数を選べば、良い計算精度が得られる
Gaussian で利用できる様々な汎関数
https://gyazo.com/151ce20fd933eaa430ef0c05764062f4
🟩 交換・相関汎関数の分類
局所スピン密度近似（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 型に変えると、ほとんどの場合、計算精度が向上
力試しテスト
密度汎関数理論（DFT）法に関する次の文章について、適切なものをすべて選んでください。
_ DFT は、HF と同様に、波動関数を試行関数とする方法であり、その基礎となる Kohn-Sham 方程式は、シュレディンガー方程式を近似したものであると言える。 x DFT は、電子同士の相互作用に由来する「電子相関」について、Hartree-Fock 法と比較して、より定量的な精度で考慮することができる計算方法である。 x 一般的に、hybrid 型の汎関数は、非 hybrid 型の汎関数に比べて、様々な物性値を計算するときに、その精度が高いと期待できる。 _ DFT の汎関数として代表的なものは B3LYP であり、この汎関数は長距離補正を含むので、分子集合体系の分子間相互作用などについても、定量的な精度で正しく評価できる。 小休憩
量子化学計算の代表的な基底関数
Gaussian で利用できる様々な基底関数
https://gyazo.com/1c1e77e0f29f800806d59bd68db3383f
基底関数とは？
任意の波動関数 $ \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
基底関数には様々な種類があり、どのような分子を解析したいのか、何を調べたいと思っているのか、どのような計算方法を用いるのか、それぞれに応じて、適切に選択する必要がある
基底関数の個数が多いほど、電子分布を柔軟に表現できるので、計算精度が高くなる
計算精度が高くなればなるほど、計算コストも高くなってしまう
計算精度と計算コストのバランスを考えて、適切に選ぶ
🟩 水素分子の場合
水素分子を表す波動関数 $ \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) レベルの計算手法と組み合わせて用いる
分極関数とは？
分子の複雑な構造を適切に記述したい場合には、各原子上に配置した基底関数に 分極関数（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）は、水素原子および水素以外の原子のどちらにも、全ての原子に分極関数が追加される
Gaussian で分極関数を追加するには
https://gyazo.com/97992beebc34a27aa6cb5f225a618e9a
分散関数とは？
負電荷を帯びた分子は、中性分子と比較して、電子分布が広がっている場合がある
このような場合、物性諸量を正しく見積もるためには、電子分布の広がりを記述するために 分散関数（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
Gaussian で分散関数を追加するには
https://gyazo.com/12865d01fe548e4315fdd25664bb973f
🟩 基底関数のベンチマーク（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）の組み合わせで、エネルギー計算を実行
力試しテスト
基底関数に関する次の文章について、適切なものをすべて選んでください。
x 最小基底系のひとつである STO-3G は、計算コストは低いが、計算精度は低いので、実験結果との定量的な比較などの実用的な解析に用いることは難しい。 x Pople 系基底関数のひとつである 6-311G は、一般に、同じ Pople 系基底関数である 3-21G よりも柔軟に電子のふるまいを記述できるため、計算精度が高い。 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) を用いた方がよい。 中休憩 15:00 〜 15:15
金属化合物の場合の注意点
金属を含む系の場合、「電子の数が多い」、特に「内殻電子が多い」ことが原因で、6-31G などでは不十分だったり、コストが高くなってしまうことがあります
このような場合、「内殻電子を近似する」ことで、適切に扱うことができる場合が多いです
具体的には...
内殻電子を「有効内殻ポテンシャル（effective core potential; ECP）」
LanL2DZ
水素原子にも分散関数を追加した方がいい場合
分子同士が水素結合した錯体を扱う場合、6-31G(d,p)などが良いかと思います
https://gyazo.com/8138f25ffada153a8f7e0a3f28c543ca
講義③ 15:15 〜 17:00
量子化学計算に基づく構造最適化
分子の安定構造の探索
各原子核に働く力（エネルギー勾配 $ \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
「カタチをちょっとずつ変化させる」を繰り返すことで、分子の立体構造を半自動的に最適化できる
実際には、多くの量子化学計算プログラムは、エネルギー勾配の方向に動かすという単純な方法（勾配降下法）ではなく、もっと効率的なアルゴリズムを使っている
Gaussian では、デフォルトでは、GEDIISというアルゴリズムを用いる
実演：Gaussian を用いた水分子の安定構造の探索
https://scrapbox.io/files/63f9a0bcd06e52001c9d4a71.mov
最適化される様子を確かめるため、HOH 結合角を 150° 程度にしたものを初期構造とした
計算例：水分子の構造最適化
https://gyazo.com/c2e381a0aff6f0bfb93bac93b78c270f
量子化学計算に基づく振動解析
分子振動とは？
水（$ \mathrm{H_2O}）の変角振動 / 対称伸縮振動 / 逆対称伸縮振動
https://gyazo.com/e4d9bb89f0f61f0533b9532d0e024287https://gyazo.com/145fab5f0795498d91045fd285f6b556https://gyazo.com/80f8eb416e3294f461be22f863b7cc58
🟧【Ｑ】化合物の反応性について、熱や光など、化合物の特性以外の要因の影響は見ることは？
「熱」については、量子化学計算では「振動解析」で得られる情報を元にして、近似的に、その影響を調べます
具体的には、分子が振動する状態を理論的に解析して、外部から「熱」が与えられたときに、分子がどのようにそのエネルギーを振動運動に変えることができるかを予測します
化学反応に伴う自由エネルギー変化など、熱力学的な性質について、近似的にですが、見積もることができます
振動解析とは？
分子固有の振動数と赤外強度を理論的に計算する方法
実験的に測定された赤外スペクトルの解釈を補助する手段として活用されている
各振動モードに対するラマン強度を計算することも可能
振動状態に基づいて、分子の熱力学的な情報（自由エネルギーなど）を計算することが可能
➡ 振動状態を解析しないと、ある温度における自由エネルギーの大きさなど、分子の熱力学的な性質についての詳しい情報が得られない
最適化した構造が「安定点」か「不安定点（遷移状態）」かを見分けることにも使える
赤外吸収スペクトルとは？
水分子の赤外吸収スペクトル（実験結果）
https://gyazo.com/246c01937cec7d08b6c54b43d41e6f3e
振動解析の計算手順
対象とする分子について、構造最適化を実行（安定構造を探索）する
構造最適化を実行したときと同じ量子化学計算レベルで、振動解析を実行する
実演：Gaussian を用いた水分子の振動解析
https://scrapbox.io/files/63f9de9406a787001b9cc00a.mov
計算例：水分子の振動解析
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{電子}
Gaussian で熱力学量を調べる
https://gyazo.com/d73533571773c87a6978e968c104c193
熱力学補正の具体例
$ \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
山本@山本研, 論文準備中
量子化学計算に基づく化学反応の解析
【応用事例】食品成分（リコピン 🍅）の trans → cis 異性化を促進したい
https://gyazo.com/84c251e60a9a9d9eeb71ddc2a807c97d
https://gyazo.com/6b1cd074ab2225a2ac3d0b11b7486e23
宮川@山本研, 分子科学討論会 (2023)
🟧【Ｑ】量子化学計算によって実験結果の予測、実験の手助けとなる計算を行うことは？
化学反応にともなう活性化エネルギー（反応障壁の大きさ）を見積もることで、「反応の速度」（反応速度定数）を理論的に見積もることができます
反応速度定数の見積もりには、「アレニウス則」や「遷移状態理論」などを用います
また、反応物と生成物のエネルギー差（反応エンタルピー）から、ある温度における「反応物と生成物の分布比」を理論的に見積もることができます
分布比は「ボルツマン分布」に従います
分子の遷移状態
ある化学反応について、反応速度を議論するときには、その反応に伴う遷移状態を調べる必要がある
たとえば、ある反応の速度定数 $ 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
【応用事例】色素材料の円偏向性を制御したい
凝集誘起発光特性をもつ三脚巴状分子の「らせん反転運動」の制御
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
遷移状態に近い（と予想される）構造から出発する
今回の場合、C-C 軸周りの回転角度を 90度くらいにしてみる
yamnor.icon 見つけようとしている遷移状態は捻れ型だと考えられるから
遷移状態探索の設定方法
Job Type
https://i.gyazo.com/c1ba79cb670331d2f93de606d041500f.png
Optimizationを選択
Optimize to a TS (Berry)
Force Constants Calculate at First Point
化学反応のエネルギーダイアグラムを描く
https://gyazo.com/4e08b3bb82bbf789a2637263df026f5d
🟩 遷移状態で振動解析すると...
値が負（虚数）となる振動数がある
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 計算ができる 小休憩
量子化学計算に基づく電子励起状態の解析
🟧【事前質問】化合物の反応性について、熱や光など、化合物の特性以外の要因の影響は見ることは？
「光」については、化合物を光で励起した後の状態（励起状態）について、量子化学計算を用いて、その影響を調べることが可能です
たとえば...
化合物が吸収・発光する光の波長 / 強度
光励起したあとの化合物の構造変化・反応過程
光の吸収とは？
光のエネルギーが分子を基底状態から励起状態に変えるエネルギーとして消費されるプロセス
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 を用いて励起状態を解析するときには、注意が必要な場合もあります。
Gaussian で TD-DFT を用いて励起状態を計算する
https://gyazo.com/d8a85658228c4b05d3ec80e566561a8c
🟩 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
量子化学計算に基づく溶媒効果の解析
🟧【Ｑ】
反応溶媒を選択する際に何か指標となるような計算ができないか？
量子化学計算で目的物の選択性を向上させるためのヒントは？
溶液中での溶解度・溶解速度や、構造の取り方（イオン化の有無 / 単一分子 o r複合体で溶解しているかなど）の予測は？
溶媒和とは
溶質分子が、静電気力や水素結合などによって結びつき取り囲むことで溶質が溶媒中に拡散する現象
https://scrapbox.io/files/653ad5ae76d858001c30094a.png
溶媒効果とは？
化学において、反応性もしくは分子の会合に対して溶媒が及ぼす影響を指す
溶媒は溶解度、安定性、反応速度に影響を及ぼすため、適切な溶媒を選択することにより、化学反応を熱力学的・速度論的に制御できる
溶媒効果の例
アセチルアセトンの互変異性に対する平衡定数 $ K_\mathrm{keto-enol} = \mathrm{\frac{[enol]}{[keto]}}
https://scrapbox.io/files/653ad6108e5379001c66ebe4.png
ある可逆的な反応過程 $ 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）
table:acetylacetone
溶媒の種類 keto-enol 平衡定数 ギブス自由エネルギー変化
気相 11.7 -6.12
水 0.23 3.67
ジクロロメタン 4.2 -3.58
エタノール 5.8 -4.39
テトラヒドロフラン 7.2 -4.93
ベンゼン 14.7 -6.70
シクロヘキサン 42 -9.32
量子化学計算で溶媒効果を考慮するには
溶媒分子を露わに考慮する（explicit solvation model）
https://scrapbox.io/files/653ad724f69f4f001b782be2.png
溶媒分子を分極する連続誘電体として扱う（implicit solvation model）
https://scrapbox.io/files/653ad7607fe76c001b737387.png
左：Debye-Onsager モデル
右：連続誘電体モデル
Gaussian で溶媒効果を考慮する
https://gyazo.com/1a93a4b7ae55da0ef53e192e745b2782
溶媒効果の具体例
エタノールについて、さまざまな溶媒に対する溶媒和エネルギー（$ \Delta E_\mathrm{sol} = E_\mathrm{sol} - E_\mathrm{gas} ）を量子化学計算に基づいて解析する
B3LYP/6-31G(d)
table:solvation effect
溶媒の種類 RHF Energy (Hartree) エネルギー差 (kcal/mol)
孤立気相 -154.075744507 ---
Water -154.081212575 -3.431
Acetonitrile -154.081069207 -3.341
Cyclohexane -154.077836155 -1.312
質疑応答 17:00 〜 17:30
質疑応答
ご質問
さいごに
共同研究・技術相談
計算化学（分子動力学・量子化学・機械学習）に関して、共同研究や技術相談をお受けしています
興味・関心を持っていただけましたら、山本までお問い合わせください
連絡先
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください