各回へのリンク 講義①
：2024-09-24講義②
：2024-10-01講義③
：2024-10-08講義④
：2024-10-15講義⑤
：2024-10-22講義⑥
：2024-10-29講義⑦
：2024-11-05講義⑧
：2024-11-12講義⑨
：2024-11-19
概要 日時： 毎週火曜日 １・２限（9:00 〜 11:00）場所： 津田沼キャンパス 653 講義室
講義資料のURL
この資料について このページでは、
山本典史
が担当する千葉工業大学の大学院講義「有機量子化学特論」の資料を限定公開 しています 講義の終了後もできる限り公開を続けたいと思っていますが、事情により、公開先を変更することがあるかも知れません
この授業について この授業では、量子化学計算 を使いはじめるときに必要となる知識と技術、量子化学計算とは何か、そしてそれを使うことで何ができるのか、その基本的な概念を説明します。 量子論や量子力学を学んでいない受講生でも、量子化学計算を実行する上で基本となる知識を身につけることができるように、分かりやすく解説します。 具体的な量子化学計算ツールとして、主に GAMESS という、無償で利用できるアプリを中心に紹介します。さらに、可視化ツールである WebMO を活用する方法も詳しく解説します。 受講生がすぐに量子化学計算に取り組めるよう、実際の入力ファイルの作成方法や計算の実行方法、結果の読み取り方なども丁寧に説明します。 この授業では具体的な例を通じて、量子化学計算の応用についても学びます。分子の安定構造 、振動状態 、励起状態 、化学反応 の解析など、実際の研究や開発で使われるような問題に取り組みながら、計算手法の実践的な使い方を身につけます。 この授業を受講する受講生には、量子化学計算をすぐに使い始めることができる WebMO のアカウント を発行します。量子化学計算のアプリをインストールする手間なく、量子化学計算を体験できます。
この講義の到達目標 🎯 量子化学計算に取り組むときに必要となる基礎的な知識と技術 について、説明できる 代表的な量子化学計算ソフトウェアである Gaussian と GAMESS の特徴と基本的な計算手順 について、説明できる 量子化学の知識や技術に基づいて、おもに有機化合物の研究開発に応用するための指針 について、説明できる 量子化学計算を活用して実験前の事前検証をおこない、研究・開発に必要となる時間的・資源的コストを効率化するための方針 について、説明できる
この授業の成績評価 中間試験・期末試験は実施しない 毎回実施する 演習課題（レポート・小テストなど） に基づいて評価をおこなう点数配分：演習課題（１３回 ）100％ manaba で「授業終了後 」に実施する場合もあるし、プリントで「授業中 」に実施することもあります
講義内容 🧑🏫 講義①： ガイダンス / 量子化学計算の概要講義②： 量子化学計算の代表的なソフトウェア講義③： 量子化学計算の基本的な手順 / 量子化学計算に基づく構造最適化量子化学計算の代表的な計算方法講義④： 基本となるハートリー・フォック法講義⑤： よく使われている密度汎関数法 講義⑥： 量子化学計算の代表的な基底関数講義⑦： 量子化学計算に基づく化学反応の解析講義⑧： 量子化学計算に基づく振動解析講義⑨： 量子化学計算に基づく電子励起状態の解析
量子化学とは？ 量子化学は、量子力学の基本原理に基づいて、分子や原子の性質や挙動を理解し説明する学問分野です 量子力学の基本原理は、シュレディンガー方程式 を基礎としています：
$$\hat{H}\Psi = E\Psi$$ ここで、$\hat{H}$ ：ハミルトニアン演算子（系の全エネルギーを表す演算子） $\Psi$ ：波動関数（系の量子状態を記述する関数） $E$ ：エネルギー固有値
量子化学計算とは？ 量子化学計算は、量子化学 を実際のコンピュータシミュレーションに適用する手法です 量子化学計算は、物質を特徴付ける微視的な性質を知りたいと思うときに、対象とする物質の 分子レベルでのふるまい を明らかにすることで、問題を解決するための手がかりを与えてくれる便利な道具です 機能性材料 を設計しようとしたり、新しい薬 を開発しようとするときに、研究・開発の現場では様々な課題に直面することもあります量子化学計算は、たとえば、物質を構成している分子の形、安定性、分子同士の相互作用の強さ、化学反応の起こりやすさなどについて、ある程度の定量的な精度で信頼できる予測 を理論的に導きだすことができます
💻 実演：量子化学計算は手軽に実行できます
量子化学計算を使ってできること 分子の安定構造の予測 振動スペクトル（IR, Raman）の予測 励起状態とUV-Visスペクトルの予測 化学反応の遷移状態と活性化エネルギーの計算 NMRケミカルシフトの計算
量子化学計算の利点は？ 実験が困難または危険な系の研究が可能 実験結果の解釈を支援 新規化合物の設計や反応経路の予測に活用可能
👋 質問：量子化学計算が「できないこと」は？ 一般に、「量子化学計算」と分類されたり、イメージされるプログラムや計算方法では、分子が数個〜十数個集まった系の計算が得意です 生体分子など、かなり大きな系でも量子化学計算ができるようになってきましたが、このような大きな系を扱おうとする場合には、立体構造をどのように準備するのか？が難しくなります多くの場合、
Amber
などの分子動力学シミュレーションなどと組み合わせて、適切な構造をサンプリングするなど、量子化学計算以外のツールも必要になります
🎈 事例： Amber を用いた分子動力学シミュレーション Mohini＠山本研 et al, PeerJ (2021)
👋 質問：化合物が固体状態でも反応予測は可能？ 化合物が固体状態の場合、もし量子化学計算で扱おうとする場合には、一部を切り出したモデル系として扱うことが多いです（「クラスターモデル」と呼ぶことがあります） 固体であることの物性（バンド構造など）が重要となる場合には、「固体物理」の分野でよく用いられている、
Quantum Espresso
などの電子状態計算プログラムをつかったほうが、解析が進むかも知れません また、動的なふるまいに興味があったり、構造の揺らぎなどを詳しく知りたいときには、
LAMMPS
などの分子動力学シミュレーションなどを利用する必要があるかも知れません
👋 質問：金属表面などへの吸着性を数値化して算出することは可能？ 金属表面への吸着を調べる場合には、上記でも少しご紹介した、クラスターモデルというアプローチで、数十個程度の金属でモデル化した表面に分子を置いて、その相互作用を調べることができます 吸着反応では、熱的な拡散過程など、動的なふるまいも重要になってきますので、量子化学計算よりも、分子動力学シミュレーションを用いることが多いかと思います分子動力学シミュレーションについて、無料で使えるものでしたら、固体表面だと LAMMPS
、多孔質材料だと RASPA
などが良いかと思います
🎈 事例： Amber を用いた分子動力学シミュレーション 井上＠山本研, 日本コンピュータ化学会（2023）
👋 質問：ある特定の分子が結晶性を有するのかどうか、計算にて考察できますか？ ある分子について、結晶構造を形成した場合の安定性を定量的に評価しようとする場合には、「
Conflex
」などを利用していただく方がスムーズかと思います分子力場モデルと呼ばれる計算手法を用いて、複数の候補構造を高速に計算することで、最も安定な結晶構造を予測することができます Gaussian で、１つの分子性結晶について詳しく量子化学計算することも可能ですCIF ファイルをそのまま読み込ませて、周期性（Periodic Boundary Condition; PBC）を考慮した計算（構造最適化など）を行うことができます ただし、一般的な量子化学計算プログラム（Gaussian や GAMESS）は「周期性」を考慮した計算が効率的ではないので、時間が掛かってしまうことが多いです
🤠 クイズ：量子化学計算を使ってできること 次に示すような「新薬を開発するケース」について、量子化学計算を用いることが適切だと思われるものをすべて選んでください。
Check
量子化学計算を体験できる Web アプリの MolCalc
量子化学エンジンは GAMESS 計算レベルは HF/STO-3G計算手法：Hartree-Fock 法 基底関数：STO-3G 開発者は Copenhagen 大学の Jan Jensen
ら 分子を SMILES 記法で指定できる
💻 実演：MolCalc を用いたエチレンの量子化学計算
SMILES 記法で分子を指定する 分子が表示されたら、Calculate Properties をクリック 計算が終わったら、調べたい性質のタブをクリック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 Frequencies3910.9 $\mathrm{cm}^{-1}$HOMO と LUMO のエネルギー差 双極子モーメントの値 Polarity and Solvation → DipoleCheck
業界標準の Gaussian 量子化学計算のデファクトスタンダード （事実上の業界標準） は Gaussian 1970年に John A. Pople
（1998年ノーベル化学賞受賞）らが最初のバージョン（Gaussian 70）を開発Pople は、
いろいろあって
、Gaussian の使用を拒否されてしまったようです 名前は、分子軌道を表現するために用いる ガウス関数（Gaussian） に由来ガウス関数：$\exp(-\alpha r^2)$ スレーター関数：$\exp(-\gamma r)$
GaussView：ユーザーインターフェース GaussView というアプリと併せて使うと、高機能な量子化学計算をお手軽に実行できる
Gaussian の特徴 計算化学者だけではなく、実験化学者にも Gaussian ユーザーが多い 分からないこと・困ったことがあるとき、適切なアドバイスをしてくれるユーザーがすぐに見つかる インターネットの検索で有益な情報をたくさん得ることができる 最新版は Gaussian 16 ライセンス方針が厳しい エネルギー計算 や 構造最適化 の収束性が良い山本は、他の量子化学計算プログラムを使うときにも、構造最適化だけは Gaussian で実行する ことがあります
Gaussian の入力ファイル 📄 例：エチレン（$\mathsf{C_2H_4}$）の構造最適化 #P HF/6-31G(d) Opt
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 を用いた量子化学計算
無料で使える GAMESS GAMESS は、Gaussian に次ぐ規模での利用実績があるG eneral A tomic and M olecular E lectronic S tructure S ystem の頭文字アメリカ版（GAMESS-US） 、イギリス版（GAMESS-UK） 、Firefly という３つの異なるバージョンがあるGAMESSと言うとき、アメリカ版（GAMESS-US） を指すことが多い どのバージョンも、1970年代後半にアメリカ・エネルギー省の NRCC（National Resource for Computational Chemistry） というプロジェクトで開発されたプログラムコードが共通の先祖 1981年オリジナル版 GAMESS が アメリカ版（GAMESS-US） と イギリス版（GAMESS-UK） に枝分かれ 1997年アメリカ版（GAMESS-US）を元にして ロシア で Firefly の開発がはじまった この３つの異なるバージョンの GAMESS は、別々のグループで開発されている各グループで独立に 、より効率的に計算できるように改良、新しい機能が追加 現在では、お互いに全く異なる性質のもの に
アメリカ 🇺🇸 版 GAMESS GAMESS-US アイオワ州立大学（アメリカ）の Mark Gordon
の 研究室
が中心に開発 学術＆商用で 無償利用可 オープンソース ではない GAMESS-US を利用した研究成果を公表する場合、開発者らによる下記の 原著論文
を引用するG. M. J. Barca, et al., J. Chem. Phys., Vol. 152, 154102 (2020) アップデートが頻繁なので、バージョン番号ではなく、更新日で管理 されている最先端で開発が進められている 様々な新しい手法 を実行できる ソースコード を無償で入手することができる
GAMESS-US の入力ファイル 📄 エチレン（$\mathsf{C_2H_4}$）の構造最適化 $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 ）について説明します
🤠 クイズ：代表的なソフトウェア 次に示す量子化学計算の代表的なソフトウェアに関する文章について、正しいと思われるものを全て選んでください。
Check
🍩 余談：Python で量子化学計算を実行してみよう Python（パイソン）は、最近注目されているプログラミング言語 たとえば Google Colab
を利用すると、お手軽に Python を試してみることができます PySCF をインストール構造最適化をするためには、geometricライブラリも必要 pip install pyscf geometric PySCF で水分子の構造最適化（HF/6-31G(d) レベル） from pyscf import gto , scf
from pyscf.geomopt.geometric_solver import optimize
mol = gto . M (
atom = "O 0.0 0.0 0.0; H 0.0 0.8 0.5; H 0.0 -0.8 0.5" ,
basis = '6-31G(d)' )
mf = scf . RHF ( mol )
opt_result = optimize ( mf )
print ( opt_result . tostring ()) O 0.00000 0.00000 -0.04868
H 0.00000 0.75468 0.52434
H -0.00000 -0.75468 0.52434
💻 実演：スパコンを用いたエタノールの量子化学計算
Gaussian や GAMESS を便利に使うための WebMO WebMO Web アプリ ブラウザーから量子化学計算を実行できる
様々な量子化学計算パッケージのインターフェースとして使えるGaussian / GAMESS / Molpro / Mopac / NWChem / ORCA / PSI4 / Quantum Espresso / Q-Chem / VASP 公式サイトの デモ版 Username: guest Password: guest 30秒 までの計算が可能タイミングが悪いと混んでいる この授業用 Username: oqc2409 Password: bn3jU4 5 分 までの計算が可能全体で 8 コア まで 分子を SMILES で指定するときには、メニューから…Lookup → Import by → SMILES
注意事項 Job Name のところに、必ず、自分の学籍番号を入力してください 自分の計算がどれか分からなくなります😢
💻 実演：WebMO を用いたエチレンの量子化学計算
🧸 実習：WebMO を用いた量子化学計算 WebMO （
授業用
）を用いて、水（H2 O）の量子化学計算（GAMESS ） を実行してください。ただし、必ず、構造最適化 （Geometry Optimization） を実行してください。計算手法には Hartree-Fock (RHF)、と基底関数には6-31G(d)を用いてください。計算が終了したら、得られた結果に基づいて、次に示す化学的性質を調べてみてください。
O-H 間の結合距離 全エネルギーの値 双極子モーメントの値 O 原子上の部分電荷の値 Check
コストを掛けずに量子化学計算をはじめる 💰 ちょっと試すGAMESS + WebMO Basic + AWS もう少し本格的にGAMESS + WebMO Pro (商用) + Sakura VPS
👋 質問：GAMESS + WebMO Basic + AWSという組み合わせを具体的に行う方法は？ 具体的なインストール方法は「
ここ
」が参考になるかと思います手順としては、クラウド（AWS や Sakura VPS）で Web サービスが実行できるように設定した後、GAMESS と WebMO Basic をインストールします 量子化学計算を実行するときには、「メモリ」を比較的多く準備する必要があります AWS の設定などは少し面倒な部分もありますので、別途、技術相談などでお話をお伺いすることもできるかと思います
量子化学計算の基本的な手順 分子の状態（電荷 とスピン多重度 ）を指定する 分子の立体構造 を指定する 分子の計算方法 と基底関数 を指定する
電荷（Charge） 計算する系全体の総電荷を指定する 例えば…１価の陽イオン（例：$\mathsf{Na}^+$） １価の陰イオン（例：$\mathsf{Cl}^-$） 一価の陽イオンと陰イオンが混在する場合（例：$\mathsf{Na}^+ + \mathsf{Cl}^-$）
スピン多重度（Multiplicity） スピン多重度は、全スピン角運動を $S$ とするとき、$2S+1$ で定義される 例えば…不対電子を持たない化合物（例：基底状態の $\mathsf{H_2}$）電子配置の例：↑↓ 全スピン角運動量： $S$ = ０ スピン多重度：$2S+1$ = １ （Singlet ) 不対電子を２個もつ化合物（例：三重項状態の $\mathsf{H_2}$）電子配置の例：↑↑ 全スピン角運動量： $S$ = １ スピン多重度：$2S+1$ = ３ （Triplet ） 不対電子を１個もつ化合物（例：$\mathsf{H_2^+}$）電子配置の例：↑ 全スピン角運動量 $S$：1/2 スピン多重度：$2S+1$ = ２ （Doublet ）
電荷とスピン多重度の例 中性分子 カチオン アニオン ラジカル 励起三重項 電荷 0 1 -1 0 0 不対電子数 0 0 0 1 2 スピン多重度 Singlet Singlet Singlet Doublet Triplet
WebMOで分子の状態を指定する
🤠 クイズ：分子の電荷とスピン多重度 次に示す分子の電荷とスピン多重度について、適切な値を答えてください。
基底状態のホルムアルデヒド $\mathsf{CH_2O}$基底状態の酸素分子 $\mathsf{O_2}$フェノールアニオン $\mathsf{C_6H_5O^-}$一酸化窒素 $\mathsf{NO}$アンモニウムカチオン $\mathsf{NH_4^+}$Check
分子の計算方法と基底関数を表記する 計算方法の種類 と 基底関数の種類 を スラッシュ （/） で挟んで表記することが多いたとえば…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で分子の計算方法と基底関数を指定する
構造最適化の基本概念 構造最適化は、分子の最も安定な立体構造を理論的に予測する重要な計算プロセスです このプロセスでは、分子のエネルギーが最小になる原子配置を探索します
キーポイント: 構造最適化の目的は、ポテンシャルエネルギー曲面上の局所最小点を見つけることです 最適化プロセスは、原子に働く力がゼロになるまで繰り返し行われます 結果として得られる構造は、実験で観測される構造と比較可能です
💻 実演：Gaussian を用いた水分子の安定構造の探索
最適化される様子を確かめるため、HOH 結合角を 150° 程度にしたものを初期構造とした
構造最適化の理論的な背景 各原子核に働く力 （エネルギー勾配 $\frac{\partial E}{\partial x}$）の方向にカタチをちょっとずつ変化（$x_\mathrm{new} = x_\mathrm{old} - \alpha \frac{\partial E}{\partial x}$）させることで、分子はより安定となる （ = エネルギーが小さくなる ）$E(x_\mathrm{new}) < E(x_\mathrm{old})$ 以下の条件を満たす原子配置を探索する$\frac{\partial E}{\partial x} = 0$
構造最適化の手順 現在の構造でエネルギーと力（エネルギー勾配）を計算 力に基づいて新しい構造を予測 新しい構造でエネルギーを再計算 収束条件を満たすまで1-3を繰り返す
安定構造を探索するシンプルな例 ある２原子分子について、原子間距離を $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.5$ である初期構造を出発して、変位の大きさを決めるパラメータ $\alpha = 0.1$ として、安定構造を探索するためのステップを進めると$r_\mathrm{new} = r_\mathrm{old} - \alpha E’(r_\mathrm{old}) = r_\mathrm{old} - 2 \alpha(r_\mathrm{old} - 1)$$= 2.5 - 2 \times 0.1 \times 1.5 = 2.2$ $E(r_\mathrm{old}=2.5) = 2.25$ かつ $E(r_\mathrm{new}=2.2) = 1.44$ なので、$E_\mathrm{new} < E_\mathrm{old}$ となる 繰り返していくと…
「カタチをちょっとずつ変化させる」を繰り返すことで、分子の立体構造を半自動的に最適化 できる 実際には、多くの量子化学計算プログラムは、エネルギー勾配の方向に動かすという単純な方法（勾配降下法 ）ではなく、もっと効率的なアルゴリズム を使っているGaussian では、デフォルトでは、GEDIIS というアルゴリズムを用いる
水分子モデルの場合
🧸 実習：水分子の構造最適化 水分子（$\mathsf{H_{2}O}$）について、
WebMO
を使って、以下の手順で構造最適化（Geometry Optimization） を実行し、結果を解析してください。新しいジョブを作成する 分子エディタで分子を描画する 計算タイプとしてGeometry Optimizationを選択する 計算方法を指定するTheory：B3LYP Basis Set：6-31G(d) ジョブを実行し、結果を調べるO-H 結合距離：0.969 Å H-O-H 結合角度：103.64 度 構造最適化で得られた結果について、実験値と比較してください。実験値：O-H 結合距離：0.958 Å H-O-H 結合角度：104.51 度 Check
構造最適化で得られた結果の解釈 最適化が完了したら、以下の点に注目して結果を解析する収束の確認 : 最適化が正常に終了したか確認構造パラメータ : 結合距離、結合角、二面角を確認エネルギー変化 : 初期構造と最終構造のエネルギー差を確認振動解析 : 最適化構造が本当に局所安定点かを確認するため、振動解析をおこなう
分子の初期配置を指定する 量子化学計算では、必ず、分子の構造 を入力する必要がある 入力する構造（初期構造 ）を出発点として、最も安定なエネルギーを持つ構造（最安定構造 ）を求める計算を「構造最適化 」（Geometry Optimization）という 構造最適化の効率は、初期構造に大きく依存する 初期構造が適切ではないと、最安定構造（global minimum）ではなく、局所安定構造 （local minimum）に収束してしまうことがある
構造最適化の注意点 計算レベルの選択 : 高精度の計算方法と大きな基底関数セットを使用すると、より正確な結果が得られますが、計算コストが増加する初期構造の重要性 : 異なる初期構造から出発すると、異なる局所最小点に到達する可能性がある
高度な構造最適化テクニック 遷移状態の最適化 : 化学反応の遷移状態を見つけるための特殊な最適化手法拘束付き最適化 : 特定の構造パラメータを固定して最適化を行う方法PES（ポテンシャルエネルギー表面）スキャン : 特定の構造パラメータを系統的に変化させて、エネルギー曲面を探索する方法
🤠 クイズ：構造最適化 構造最適化に関する以下の記述のうち、正しいものをすべて選んでください。
Check
多電子系のシュレディンガー方程式 多電子系の場合、シュレディンガー方程式は以下のように展開されます：
$$\left[ -\frac{1}{2} \sum_i^n \nabla_i^2 + \sum_i^n v(r_i) + \sum_{i}^n \sum_{j > i}^n \frac{1}{\left| r_{i} - r_{j} \right|} \right] \Psi = E \Psi$$ ここで、第1項：電子の運動エネルギー 第2項：電子-核間の引力 第3項：電子-電子間の反発 しかし、この方程式を厳密に解くことは、水素原子のような最も単純な系を除いて、実質的に不可能です そこで、多電子系のシュレディンガー方程式を実用的に解くために、量子化学計算の様々な近似手法が開発されました
量子化学計算の計算手法 電子相関を考慮しない計算方法 電子相関を考慮する計算方法密度汎関数法 Møller-Plesset 摂動法（例：MP2, MP4） 結合クラスター法（例：CCSD, CCSD(T)） 配置間相互作用（例：CISD, Full CI）法
💻 実演：Gaussian で利用できる様々な計算方法
シュレディンガー方程式 $$
\left[ -\frac{1}{2} \sum_i^n \nabla_i^2 + \sum_i^n v(r_i) + \sum_{i}^n \sum_{j > i}^n \frac{1}{\left| r_{i} - r_{j} \right|} \right] \Psi = E \Psi
$$
２項目の $v(r_i)$ は、$i$ 番目の電子に外部（原子核）からはたらくポテンシャル
$$
v(r_i) = -\sum_a^{N_\mathbf{nuclei}} \frac{Z_a}{|R_a - r_i|}
$$ ３項目の $\frac{1}{|r_{i}-r_{j}|}$ は、$i$ 番目と $j$ 番目の「電子間にはたらく相互作用（電子相関） 」を表す「ある電子は、他の電子の位置に依存しながら運動している 」という描像（多電子描像 ）厳密に取り扱うことが難しいため、何らかの「近似 」が必要
Hartree-Fock (HF) 方程式 量子化学計算の出発点となる最も基本的な近似 「ある電子は、他の電子の作る平均化された場の中で（他の電子と）独立に運動している 」という描像（一電子近似 ）を考える
$$
\left[ -\frac{1}{2} \nabla_i^2 + v(r_i) + V_i^{\mathrm HF} \right] \psi_i = \epsilon_i \psi_i
$$
３項目の $V_i^{\mathrm HF}$ は、$i$ 番目の電子に対して他の電子（$j \neq i$）が作る「平均場 」に由来する反発エネルギー
🟩 ステップ 1️⃣ 対象とする分子系のハミルトニアン（核座標の配置 $R_\alpha$ 、核電荷 $Z_\alpha$、電子数 $n$ などの情報）、初期の参照軌道 $\Psi^{0}$ を設定する
対象とする分子系のハミルトニアン：
$$
\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}(r_1)\;\psi^{0}(r_2)
$$もちろん、ここにある $\psi^0(r_1)$ や $\psi^0(r_2)$ は、「これから求めたいもの」なので、その中身が分からない 拡張 Huckel 近似や半経験法など、精度の粗い方法をつかって、とりあえず求めておく
🟩 ステップ 2️⃣ 「電子２がつくる平均場」中にある電子１について、HF 方程式を解く
$$
\left[-\frac{1}{2}\nabla_1^2 - \frac{2}{r_1} + \int d r_2 \; |\psi^0(r_2)|^2 \, \frac{1}{r_{12}}\right]\psi(r_1)=\epsilon_1 \psi(r_1)
$$
これを解くと電子１の波動関数 $\psi(r_1)$ とエネルギー $\epsilon_1$ が求まるけれど、ステップ①で「とりあえず」ということで決めた $\psi^0(r_2)$ をつかった計算なので、精度的にはまだまだ不十分
🟩 ステップ 3️⃣ 「電子１がつくる平均場」中にある電子２について、HF 方程式を解く
$$
\left[-\frac{1}{2}\nabla_2^2 - \frac{2}{r_2} + \int d r_1 \; |\psi^0(r_1)|^2 \, \frac{1}{r_{12}}\right]\psi(r_2)=\epsilon_2 \psi(r_2)
$$
これを解くと電子２の波動関数 $\psi(r_2)$ とエネルギー $\epsilon_2$ が求まるけれど、、ステップ①で「とりあえず」ということで決めた $\psi^0(r_1)$ をつかった計算なので、得られる答えも、精度的にはまだまだ不十分
🟩 ステップ 4️⃣ 軌道エネルギー $\epsilon_i$ と１電子軌道 $\psi(r_i)$ から、全電子エネルギー $E$ を求める
ステップ②の続きでは、アップデートした $\psi^0(r_2)$ をつかった計算なので、得られる答えは、１回目よりも改善される（エネルギーが下がる） 得られた全エネルギー $E_\mathrm{new}$ と前回のもの $E_\mathrm{old}$ とのエネルギー差 $\Delta E$ を比べたときに…$\Delta E$ が十分に小さい $\Delta E$ が大きい➡ アップデートした軌道 $\Psi_\mathrm{new}$ を参照 $\Psi^0$ に使って、ステップ②〜④を再度実行
水素分子の場合
Hartree-Fock 法の特徴 分子の立体構造や化学的性質について、多くの場合、十分な精度で予測 できる分子が最適構造にあるとき、その分子の全エネルギーの値を99%程度の精度 で見積もることができる しかし、化学反応を議論するためには、99%程度という精度でも不十分 であり、より高精度な予測が必要となる化学反応に伴う活性化エネルギー $E_a$ を議論しようとする場合、$E_a$ は 20 kJ/mol 程度（5 kcal/mol 程度） になることが多いので、この程度以下の誤差となることが望ましい 電子相関 を考慮するように HF 法を拡張し、HF 法以上の計算精度を持つ手法のことを「post Hartree-Fock 法 」と呼ぶpost HF 法では、HF 法では無視してしまった「電子相関 」の効果について、様々な工夫で考慮する
🤠 クイズ：Hartree-Fock 法 Hartree-Fock (HF) 法に関する次の文章について、適切なものをすべて選んでください。
Check
電子相関とは？ Hartree-Fock 法は、電子同士の相互作用（電子相関 ）を平均化して扱うという近似（平均場近似 ）をしている
電子相関の大きさは「電子相関エネルギー 」（$\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）とは？ 近年、研究開発の最前線で大活躍している計算手法 電子密度を試行関数 とする多電子状態問題の近似解法高田 (1999)
電子相関の寄与について、 Hartree-Fock 法とほぼ同程度の計算コストで考慮 することができるHartree-Fock 法 などは波動関数を試行関数 とする方法であるが、DFT は電子密度を試行関数 とする方法であり、その基盤となる考え方が大きく異なっている HF 法の基礎は HF 方程式 だが、DFTを実践する上で鍵になるのは コーン・シャム（Kohn-Sham; KS）方程式
Hartree-Fock 方程式 $$
\left[ -\frac{1}{2} \nabla_i^2 + v(r_i) + V_i^\mathrm{HF} \right] \psi_i = \epsilon_i \psi_i
$$
３項目の $V^\mathrm{HF}$ は「ある電子-他の電子間の「平均」反発エネルギー 」を表す 「ある電子は、他の電子の作る平均化された場の中で（他の電子と）独立に運動している」という描像
Kohn-Sham 方程式 $$
\left[ -\frac{1}{2} \nabla_i^2 + v(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}$：非古典的な相互作用であり、（交換相関）汎関数と呼ばれる このアイデアを大まかに説明すると…外部有効ポテンシャル $V^\mathrm{eff}$ のもとにある $n$ 個の相互作用がない仮想の電子系 （Kohn-Sham 系）を考える ただし $V^\mathrm{eff}$ は、相互作用を考慮したリアルな系で求められるものと等しい（仮想的な分子系） とする
Hartree-Fock 法と DFT 法（Kohn-Sham 法）の違い Kohn-Sham 方程式と Hartree-Fock 方程式の違いは、３項目が「平均場 $V^\mathrm{HF}$」か「外場有効ポテンシャル $V^\mathrm{eff}$」かのみつまり、Kohn-Sham 方程式は、Hartree-Fock 方程式と同程度の時間で計算できる Hartree-Fock 方程式中の「平均場 $V^\mathrm{HF}$」は「既知の関数 」であるが、Kohn-Sham 方程式中の「交換相関汎関数 $V_\mathrm{xc}$」は厳密な形が不明な「未知の関数 」である$V_\mathrm{xc}$ の厳密な形が分からないので、人為的に仮定された汎関数 を使う この 「人為的に仮定された」汎関数の質 が、DFT の精度を決定する
密度汎関数法で用いる「汎関数」 DFT の計算精度と信頼性は、電子密度と系のエネルギーを関係付ける汎関数 の種類に依存する 適切な汎関数を選べば、良い計算精度が得られる
Gaussian で利用できる様々な汎関数
🟩 交換・相関汎関数の分類 局所スピン密度近似（LSDA）型 一般化勾配近似（GGA）型 LSDA を電子密度勾配 $\nabla \rho$ で補正した汎関数 meta GGA 型 GGA を運動エネルギー密度 $\tau$で補正した汎関数 hybrid 型 Hartree-Fock 交換積分 $E_\mathrm{x}^\mathrm{HF}$ を一定量混合した汎関数
基本的な GGA 汎関数 BLYP （Becke-Lee-Yang-Parr） PBE（Perdew-Burke-Ernzerhof）
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 など）
🧸 実習：水分子の計算 水分子（$\mathsf{H_2O}$）について、
WebMO
を使って、下記で指定される計算手法で構造最適化（Geometry Optimization） を実行してください。Theory（※下記からひとつ好きなものを選んで）:DFT BLYPDFT B3LYPDFT PBEDFT PBE0 Basis Set: 計算結果から、O-H 結合距離、および、H-O-H結合角度の値をそれぞれ読み取り、次に示す実験値と比較してください。実験値：O-H 結合距離：0.958 Å H-O-H 結合角度：104.51 度 上記が終わったら、計算手法（Theory / Basis Set）を変えて、計算結果を比較してみてください。
計算例：水分子の構造最適化
hybrid 交換・相関汎関数 📊 hybrid 型は、LSDA / GGA / mGGA などの交換汎関数に対して、一定の割合で Hartree-Fock 交換積分 $E_\mathrm{x}^\mathrm{HF}$ を混合した汎関数です。 非 hybrid 型（local 型）と hybrid 型の汎関数ペア（例：BLYP vs B3LYP）について、計算される物性値の平均二乗偏差（RMSD）を比較してみると…NCD：分子間相互作用エネルギーID：異性化エネルギーBH：反応障壁の高さ Mardirossian & Head-Gordon (2017)
$E_\mathrm{x}^\mathrm{HF}$ を適度に hybrid することで計算精度が向上
分散力補正 📊 標準的な DFT の弱点のひとつは、分散力（van der Waals 相互作用） の記述ができないこのために、分子集合体における分子間相互作用を過小評価する傾向にある Grimme は、分子力場などで用いられるように経験的なパラメータを用いて、各原子対に対する引力的なエネルギー項を付加する分散力補正 を提案 Grimme の分散力補正のあり・なし（例：BLYP と BLYP-D2, BLYP-D3）について、計算される物性値の平均二乗偏差（RMSD）を比較してみると…NCED：分子間相互作用エネルギーIE：異性化エネルギー Mardirossian & Head-Gordon (2017)
Grimme の分散力補正を加えることで、コストを追加することなく、計算精度を向上させることができる
交換・相関汎関数のベンチマーク 📊 200 個の汎関数について、ベンチマーク計算をしてみると…汎関数名の右にある数字は、計算精度のランキング 左にあるラベルは、評価した物性値：NCED / NCED / NCD：分子間相互作用エネルギーIE / ID：異性化エネルギーTCE / TCD：熱化学量BH：反応障壁の高さ Mardirossian & Head-Gordon (2017)
GGA から meta GGA に変えても、多くの場合、計算精度はそれほど良くはならない 非 hybrid 型（local 型）から hybrid 型に変えると、ほとんどの場合、計算精度が向上
🤠 クイズ：密度汎関数理論法 密度汎関数理論（DFT）法に関する次の文章について、適切なものをすべて選んでください。
Check
Gaussian で利用できる様々な基底関数
基底関数とは？ 任意の波動関数 $\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}$ の値だけを変える 。
基底関数の選び方 基底関数には様々な種類があり、どのような分子 を解析したいのか、何を調べたい と思っているのか、どのような計算方法 を用いるのか、それぞれに応じて、適切に選択する必要がある基底関数の個数が多いほど、電子分布を柔軟に表現できるので、計算精度 が高くなる 計算精度 が高くなればなるほど、計算コスト も高くなってしまう計算精度 と計算コスト のバランスを考えて、適切に選ぶ
🟩 水素分子の場合 水素分子を表す波動関数 $\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)
$$
したがって、水素原子の 1s 軌道を参考にして、$\exp(-\zeta r)$ 型の関数を基底関数として使っても良さそうだ $$
\phi_a = \exp(-\zeta r)
$$
このような関数をスレーター型関数 と呼ぶ 係数 $\zeta$ は、関数の「広がりの程度」を表す変数 水素原子の厳密な解析解（= 電子の分布を適切に表すことができる）なので、計算精度は高い
ただし、数値計算の計算効率は悪い 図にあるように、原子の中心（$r = 0$）で値が不連続になるから、数値的に積分することが難しい 計算効率が悪いので、ふつうは、量子化学計算では使わない
🟩 ガウス型関数 量子化学計算では、スレーター型関数 $\exp(-\zeta r)$ の代わりに、$\exp(-\alpha r^2)$ 型の関数を基底関数として使う $$
\phi_a = \exp(-\alpha r^2)
$$
係数 $\alpha$ は、関数の「広がりの程度 」を表す変数量子化学計算の前に、あらかじめ設定されている定数（計算の途中で変化させない） ガウス型関数は、原子軌道に対する近似解 なので、計算精度は高くはない ただし、数値計算の計算効率が良い 上の図にあるように、原子の中心（$r = 0$）でも値が連続になるから、数値積分が簡単 計算精度の悪さは、拡がりの程度を変えたガウス型の基底関数をいくつか線形結合させる ことでカバーするいくつかのガウス型関数を線形結合させたものを# 短縮型ガウス基底とよぶ たとえば… $$
\phi_\mathsf{blue} = c_\mathsf{cyan} \, \chi_\mathsf{cyan} + c_\mathsf{pink} \, \chi_\mathsf{pink} + c_\mathsf{green} \, \chi_\mathsf{green}
$$
$$
= \left(1.33 \; e^{-0.59 \; r^2_\mathsf{cyan}}\right) + \left(-0.57 \; e^{-1.28 \; r^2_\mathsf{pink}}\right) + \left(0.04 \; e^{+0.06 \; r^2_\mathsf{green}}\right)
$$
青い線で示す 短縮型ガウス関数 は、赤い線で示すスレーター型関数のかたちを、まあまあ良く再現している
最小基底系（STO-nG基底系） 例STO-3G３種類 のガウス型 関数（3G）を足し合わせて、スレーター型 軌道（STO）に似せたものSTO-4GSTO-5G 計算精度の比較 計算コストは低いが、計算精度も低いので、実験結果と比較するなどの実用的な解析に用いることはあまりない
Pople 系 例 計算精度の比較 John A. Pople
（1998年ノーベル化学賞受賞）らが提案した基底関数系最小基底系に比べると、より柔軟に電子のふるまいを記述できる ため、計算精度が改善される原子価軌道に対して、拡がりの異なる複数の関数を組み合わせて表現する 一般的な有機化合物に対しては、基盤となる 6-31G に 分極関数 を加えた 6-31G(d) や 6-31G(d,p) などが、計算効率が良い・実用的な精度をもつ基底関数としてよく用いられている
分割価電子（split valence）基底関数 内殻軌道には１つ、価電子軌道には複数の短縮ガウス型関数を用いた基底関数系 Pople 系の場合
6-31G の場合 １つの占有軌道に対して…内殻軌道６ 個のガウス型を足し合わせた（短縮した）１種類 の関数で記述 価電子軌道３ 個のガウス型を短縮したもの＋別の１ 個のガウス型を組み合わせて、計２種類 の関数で記述２種類の関数で原子価軌道を表すものを double zeta (DZ) 型とよぶ
6-311G の場合 １つの占有軌道に対して…内殻軌道６ 個のガウス型を足し合わせた（短縮した）１種類 の関数で記述 価電子軌道３ 個のガウス型を短縮したもの＋別の１ 個のガウス型＋別の１ 個のガウス型を組み合わせて、計３種類 の関数で記述３種類の関数で価電子軌道を表すものをtriple zeta (TZ) 型とよぶ
具体例
高精度計算に適した基底関数（Correlation-Consistent 基底系） 電子相関を取り込んだ高精度計算に適した基底関数として Thom H. Dunning, Jr.
らの cc-pVnZ 基底関数（n = D, T, Q, 5, 6）がある 計算精度の比較cc-pVDZ < cc-pVTZ < cc-pVQZ 応用研究の中では、double zeta 精度の cc-pVDZ が頻繁に用いられるdouble zeta 精度では、電子相関効果の取り込みが不十分である場合もある化学反応などにおいて 10 kcal/mol 程度のエネルギー差を定量的に議論したいときには、cc-pVTZ 程度（triple zeta 精度）の基底関数を CCSD(T) レベルの計算手法と組み合わせて用いる
分極関数とは？ 分子の複雑な構造を適切に記述したい場合には、各原子上に配置した基底関数に 分極関数（polarization function） を追加することで、電子分布が柔軟に変形する自由度を与えることができる
分極関数（d 型関数） p 軌道をもつ水素以外の原子 に対しては、d 型の関数を追加する
分極関数（p 型関数） 球対称な s 軌道のみを持つ水素原子 に対しては、p 型の関数を追加する
表記方法 / 指定方法 🔤 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 で分極関数を追加するには
分散関数とは？ 負電荷を帯びた分子 は、中性分子と比較して、電子分布が広がっている 場合があるこのような場合、物性諸量を正しく見積もるためには、電子分布の広がりを記述する ために 分散関数（diffuse function） を加えることができる
表記方法 / 指定方法 🔤 Pople 系に対して、分子内にある水素以外の原子 に分散関数 を追加する場合3-21+G(d) or 3-21+G*6-31+G(d) or 6-31+G* Pople 系に対して、水素以外の原子 に分散関数 を追加した上で、水素原子にも分散関数を追加する場合 cc-pVnZ 基底関数（n = D, T, Q, 5, 6）に対して、水素原子および水素以外の原子のどちらにも、全ての原子 に分散関数 を追加する場合
Gaussian で分散関数を追加するには
🤠 クイズ：基底関数 基底関数に関する次の文章について、適切なものをすべて選んでください。
Check
オススメの計算手法と基底関数は？ 一般的な有機化合物の基底状態の物性を調べようとする場合、まずは B3LYP/6-31G(d) や B3LYP/6-31G(d,p) を使うことが多いです原子間の結合距離 を定量的に議論するには、最低でも double zeta 精度 （6-31G）が必要 原子間の結合角度 を定量的に議論するには、分極関数 （6-31G(d)）が必要
🟩 使いたい基底関数がプログラムに実装されていない場合は？ %chk=methane.chk
#P HF/Gen
Methane
0 1
C 0.000000 0.000000 0.000000
H 0.000000 0.000000 1.089000
H 1.026719 0.000000 -0.363000
H -0.513360 0.889165 -0.363000
H -0.513360 -0.889165 -0.363000
H 0
S 2 1.00
0.5447178000D+01 0.1562849787D+00
0.8245472400D+00 0.9046908767D+00
S 1 1.00
0.1831915800D+00 1.0000000
****
C 0
S 6 1.00
0.3047520000D+04 0.1825880123D-02
0.4564240000D+03 0.1405660094D-01
0.1036530000D+03 0.6875700462D-01
0.2922580000D+02 0.2304220155D+00
0.9348630000D+01 0.4684630315D+00
0.3189040000D+01 0.3627800244D+00
SP 2 1.00
0.3664980000D+01 -0.3958951621D+00 0.2364599466D+00
0.7705450000D+00 0.1215834356D+01 0.8606188057D+00
SP 1 1.00
0.1958570000D+00 0.1000000000D+01 0.1000000000D+01
****
キーポイント 量子化学計算を用いることで、実験では直接観測が困難な反応中間体や遷移状態の性質を理解することができます 以下のような化学反応に関する重要な情報を得ることができます：反応物と生成物のエネルギー差（反応エンタルピー） 反応経路上の遷移状態の構造とエネルギー
🎈 事例： 食品成分（リコピン 🍅）の異性化を促進したい 宮川@山本研, 分子科学討論会 (2023)
💻 実演：Gaussian を用いた Diels-Alder 反応の解析
👋 質問：量子化学計算によって実験結果の予測、実験の手助けとなる計算を行うことは？ 化学反応にともなう活性化エネルギー （反応障壁の大きさ ）を見積もることで、「反応の速度 」（反応速度定数 ）を理論的に見積もることができます反応速度定数の見積もりには、「アレニウス則 」や「遷移状態理論 」などを用います また、反応物と生成物のエネルギー差（反応エンタルピー ）から、ある温度における「反応物と生成物の分布比 」を理論的に見積もることができます
分子の遷移状態 ある化学反応について、反応速度を議論するときには、その反応に伴う 遷移状態（Transition State; TS） を調べる必要がある
遷移状態の基本概念 定義 反応経路上で最もエネルギーが高い点 反応物から生成物へ変化する際の中間的な構造 エネルギー的特徴 ポテンシャルエネルギー曲面上のサドルポイント 一次の鞍点：1つの方向でエネルギー極大、他の全ての方向で極小 構造的特徴 反応物と生成物の中間的な構造 形成中/開裂中の結合が通常の結合長と異なる 正反応と逆反応の遷移状態は同一 Hammond の仮説 遷移状態の構造は、エネルギー的に近い状態（反応物 or 生成物）に類似 実験的検出の難しさ 寿命が極めて短い（フェムト秒オーダー） 直接観測は困難で、計算化学が重要な役割を果たす 振動解析での特徴 1つの虚振動（負の振動数）を持つ この虚振動が反応座標に対応
反応速度を予測するには？ たとえば、ある反応の速度定数 $k$ について、始状態と遷移状態とのエネルギー差 、つまり、活性化エネルギー $\Delta G^{\ddagger}$ が分かっていれば、Eyring-Polanyi 式 から…$\chi$：遷移状態において、反応座標に沿った分子振動が一往復したときに反応する確率（透過係数 ） $h$：
プランク定数
遷移状態に基づく解析の具体例 $\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}$ 重水素 （$\mathsf{D}$） に置換した場合$\Delta G^{\ddagger}$= 18.86 kcal/mol $k(298\mathsf{K})$ = 0.0928 $\mathsf{s}^{-1}$
🎈 事例： 色素材料の円偏向性を制御したい 大場@山本研, 分子科学討論会 (2023)
遷移状態の探索 多くの量子化学計算プログラムには、遷移状態を探索する機能 が備わっている単純に「構造最適化の逆」（$x_\mathrm{new} = x_\mathrm{old} + \alpha \frac{\partial E}{\partial Q}$）をやっても、うまくいかない
エネルギーの二次微分（$\frac{\partial^2 E}{\partial Q^2}$）を要素にもつ Hessian 行列 の固有値（力の定数 or 曲率：$\chi$）を求めて…Hessian 行列の負の固有値を持つ固有ベクトルの方向 に少しずつ動かすという方法がある→ Eigenvector Following 法 調和近似（$E = \kappa Q^2$）からの歪みが大きな方向に少しずつ動かすという方法がある Hessian 行列を計算するためのコスト がまあまあ大きい
👋 質問：反応の場などの前提条件は、事前に設定したうえで遷移状態を探索するのでしょうか？ あらかじめ、「遷移状態に近いのはどのような構造か？」について、大まかな知識や前提が必要となります
👋 質問：遷移状態を見つけるのが難しい場合には？ Nudged Elastic Band (NEB) 法や string 法
などを用いて、最小エネルギー経路（Minimum Energy Path）を探索する方法があります
💻 実演：Gaussian を用いたブタジエンの cis/trans 異性化反応の解析
遷移状態に近い（と予想される）構造から出発する今回の場合、C-C 軸周りの回転角度を 90度くらいにしてみる見つけようとしている遷移状態は捻れ型だと考えられるから 遷移状態探索の設定方法
Optimization を選択Optimize to a TS (Berry) Force Constants Calculate at First Point
🤠 クイズ：化学反応の解析 ブタジエンの cis/trans 異性化反応について、量子化学計算で得られた結果は以下の通りです：trans 体のエネルギー $\epsilon_\mathrm{trans} =$ -155.992139 Hartree cis 体のエネルギー $\epsilon_\mathrm{cis} =$ -155.986486 Hartree 遷移状態のエネルギー $\epsilon_\mathrm{TS} =$ -155.980090 Hartree 量子化学計算で得られた結果に基づいて、以下の問いに答えてください。ただし、気体定数の値は $R = 1.987206 \ \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.987206 \ \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.987206 \ \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.987206 \ \mathrm{cal \ K^{-1} \ mol^{-1}} \ \times \ 300 \ \mathrm{K}}\right)$ = 1.192 $\times \ 10^{11} \ \mathrm{sec^{-1}}$ Check
🤠 クイズ：化学反応の解析 量子化学計算に基づく化学反応の解析に関する各問について、正しい選択肢を選んでください。
問１） 遷移状態の特徴として正しいものは？
問２） 活性化エネルギーの定義として正しいものは？
問３） 遷移状態の振動解析で期待される結果は？
問４） Hammond の仮説によると、遷移状態の構造は？
問５） 反応エネルギーが負の値の場合、その反応は？
分子振動とは？ 水（$\mathrm{H_2O}$）の変角振動 / 対称伸縮振動 / 逆対称伸縮振動
分子振動の自由度 分子振動の自由度（振動モードの数） $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$
🤠 クイズ：振動自由度の数 水（H2 O）、二酸化炭素（CO2 ）、メタノール（CH3 OH）の振動自由度の数はいくつですか？
水：3 × 3 - 6 = 3 二酸化炭素：3 × 3 - 5 = 4 メタノール：3 × 6 - 6 = 12 Check
👋 質問：化合物の反応性について、熱や光など、化合物の特性以外の要因の影響は見ることは？ 「熱 」については、量子化学計算では「振動解析 」で得られる情報を元にして、近似的に 、その影響を調べます具体的には、分子が振動する状態を理論的に解析して、外部から「熱」が与えられたときに、分子がどのようにそのエネルギーを振動運動に変えることができるかを予測します 化学反応に伴う自由エネルギー変化など、熱力学的な性質について、近似的にですが、見積もることができます
振動解析とは？ 分子固有の振動数と赤外強度 を理論的に計算する方法 実験的に測定された赤外スペクトル の解釈を補助する手段として活用されている 各振動モードに対するラマン強度 を計算することも可能 振動状態に基づいて、分子の熱力学的な情報 （自由エネルギーなど）を計算することが可能➡ 振動状態を解析しないと、ある温度における自由エネルギーの大きさなど、分子の熱力学的な性質についての詳しい情報が得られない 最適化した構造が「安定点」か「不安定点（遷移状態）」かを見分けることにも使える
振動解析の計算手順 対象とする分子について、構造最適化を実行 （安定構造を探索）する 構造最適化を実行したときと同じ量子化学計算レベル で、振動解析を実行する
💻 実演：Gaussian を用いた水分子の振動解析
🧸 実習：水分子の振動解析 水分子（H2 O）について、
WebMO
を使って、下記で指定される計算手法で振動解析（Vibrational Frequencies） を実行してください。その際、必ず、振動解析の前に「構造最適化（Geometry Optimization） 」をおこなってください。Theory（※下記から好きなものを選んで）:DFT BLYPDFT B3LYPDFT PBEDFT PBE0 Basis Set（※下記から好きなものを選んで）: 計算結果から、変角振動、対称伸縮振動、逆対称伸縮振動の値をそれぞれ読み取り、次に示す実験値との誤差を計算してください。実験値：変角振動：1595 cm-1 対称伸縮振動：3657 cm-1 逆対称伸縮振動：3756 cm-1 上記が終わったら、計算手法（Theory / Basis Set）を変えて、計算結果、および、計算に掛かった時間を比較してみてください。
計算例：水分子の振動解析
遷移状態で振動解析すると…
安定構造では… エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロ である
$$
\frac{\partial E}{\partial Q} = 0
$$局所 安定状態（下に凸の放物線） であり、エネルギーの座標に関する二次微分の値が常に正 である
$$
\frac{\partial^2 E}{\partial Q^2} > 0
$$
遷移状態では… エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロ である
$$
\frac{\partial E}{\partial Q} = 0
$$特定の変位方向に対して局所 不安定状態（上に凸の放物線） であり、その変位方向にはエネルギーの座標に関する二次微分の値が負 となる
$$
\frac{\partial^2 E}{\partial Q^2} < 0
$$ 振動解析すると、１つの振動モードの振動数のみが負の値になる
分子振動の非調和性
MP2摂動法を用いることで、実験値にかなり近くなる 高波数の振動モード（対称伸縮＆逆対称伸縮）では、誤差が大きい分子振動を 調和振動子
だと見なす近似（調和近似 ）による誤差（非調和性）の影響だと考えられる
スケーリング因子 非調和性の問題を手軽に解決する場合には、振動数に適当な数字（スケーリング因子 ）を掛けて補正する
🟩 非調和振動数解析 分子振動をより定量的精度で解析したい場合には、調和近似を越えた取り扱いが必要 GAMESS では、非調和性を考慮した振動解析として、Vibrational SCF (VSCF) 法 を利用することができる VSCF法で計算した場合
振動解析に基づいた熱力学量の見積もり 分子の振動状態に基づいて、電子状態のエネルギーに対して、 熱力学量を「近似的に」算出するための補正 をおこなうことができる振動数が低い振動モードほど、熱力学補正に対する寄与が大きい 「近似」なので、大きく構造が揺らぐ柔軟な分子については、誤差が大きくなる場合がある 求めることができるのは…ゼロ点エネルギー（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 で熱力学量を調べる
GAMESS で熱力学量を調べる -------------------------------
THERMOCHEMISTRY AT T= 298.15 K
-------------------------------
（省略）
E H G CV CP S
KCAL/MOL KCAL/MOL KCAL/MOL CAL/MOL-K CAL/MOL-K CAL/MOL-K
ELEC. 0.000 0.000 0.000 0.000 0.000 0.000
TRANS. 0.889 1.481 -8.837 2.981 4.968 34.608
ROT. 0.889 0.889 -2.615 2.981 2.981 11.753
VIB. 14.419 14.419 14.418 0.023 0.023 0.003
TOTAL 16.196 16.789 2.965 5.985 7.972 46.364
VIB. THERMAL CORRECTION E(T)-E(0) = H(T)-H(0) = 0.776 CAL/MOL
熱力学補正の具体例 $\mathsf{C_2H_5 + H_2 \rightarrow C_2H_6 + H}$ Wang and Zhao, PCCP (2011)
エネルギー変化（熱力学補正をしない場合）$\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
🤠 クイズ：振動解析 量子化学計算に基づく振動解析に関する各問いについて、正しい選択肢を選んでください。
問１） メタン（CH4 ）の振動モードの数は？
問２） 振動解析で得られる負の振動数は通常何を示唆しますか？
問３） 振動解析結果から直接計算できる熱力学量はどれですか？
問４） 同位体置換（例：H2 O → D2 O）が振動スペクトルに与える主な影響は？
問５） スケーリング因子を振動数に適用する主な目的は何ですか？
問６） 赤外不活性だがラマン活性である振動モードの特徴は何ですか？
👋 質問：化合物の反応性について、熱や光など、化合物の特性以外の要因の影響は見ることは？ 「光 」については、化合物を光で励起した後の状態（励起状態）について、量子化学計算を用いて、その影響を調べることが可能です たとえば…化合物が吸収・発光する光の波長 / 強度 光励起したあとの化合物の構造変化・反応過程
光の吸収とは？ 光のエネルギーが分子を基底状態から励起状態に変えるエネルギー として消費されるプロセス
$E_{励起} - E_{基底} = \Delta E = h c / \lambda_{abs}$ に相当する光のエネルギーを分子が吸収し、基底 → 励起状態への電子遷移 が起こる
分子の吸光・発光と基底・励起状態
吸光 基底（$\mathrm{S_0}$）状態の最適構造 → 励起（$\mathrm{S_n}$）状態つまり、量子化学計算に取り組むときに、分子の立体構造は、基底状態で最適化する必要がある 発光 励起（$\mathrm{S_n}$）状態の最適構造 → 基底（$\mathrm{S_0}$）状態つまり、量子化学計算に取り組むときに、分子の立体構造は、励起状態で最適化する必要がある
🎈 事例： 高効率な有機ＥＬを設計したい 分子軌道のカタチ を調べることでロジカルな設計が可能となった➡ HOMO と LUMO の重なりを少なくすることで、$\mathsf{T_1}$ → $\mathsf{S_1}$ の項間交差（ISC）を促進する 安達＠九大, et al., Nature Comm., Vol. 6, p. 8476 (2015)
🎈 事例： 環境応答性をもつ色素材料を設計したい 古川@山本研, 分子科学討論会 (2023)
励起状態の計算方法 TD-DFT 法 時間依存密度汎関数理論（Time-Dependent DFT）法。励起状態の計算方法として、最も普及している CIS 法 SAC-CI 法 Symmetry Adapted Cluster/Configuration Interaction (SAC-CI) 法 CASSCF 法 完全活性空間 (Complete Active Space) SCF 法
DFT を用いて励起状態を解析できる？ DFT は、基底状態の計算に関しては、HF 法と同等のコストで、post HF 精度 の結果が得られる可能性があります。 DFT は、原子や分子の励起状態 に関しても、時間依存（TD）の問題として扱うように拡張（TD-DFT ）できます。TD-DFTは、化学発光 / 光触媒 / 太陽電池 / 光合成など、さまざまな分野で活用されています。 ただし、TD-DFT を用いて励起状態を解析するときには、注意が必要な場合もあります。
Gaussian で 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 方程式を解くと…
$$
\begin{bmatrix} A & B \\ B^* & A^* \end{bmatrix}
\begin{bmatrix} X \\ Y \end{bmatrix}
= \omega \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix}
\begin{bmatrix} X \\ Y \end{bmatrix}
$$
$\pm \omega$ の値が、励起 ＆ 脱励起のエネルギーに対応します
行列 $\mathbf{B}$ を無視して上記の方程式を解く扱いのことを、Tamm-Dancoff 近似（TDA） といいます：
$\mathbf{AX} = \omega \mathbf{X}$ 遷移モーメントの定量性が悪くなる可能性があります 一方で、励起エネルギーの計算精度が改善されることがあります
TD-DFT のベンチマーク（TDA など） GINV：meta GGA 汎関数 $E_\mathrm{xc}^{mGGA}[\rho,\nabla\rho,\tau]$ の運動エネルギー密度 $\tau(\mathbf{r},t)$ のゲージ依存性を補正した結果Liang, et al. (2022) 三重項状態の計算では、TDA によって計算精度が改善しています。 TDA はゲージ不変性を破ってしまうので、GINV と同時に採用することはできません。GINV による改善はわずかなので、通常は、TDA の方をオススメします。
TD-DFT のベンチマーク（汎関数依存性） 📊 Liang, et al. (2022)
領域分割 hybrid 型（長距離補正）を用いることで、計算精度が改善しています。 蛍光（fluorescence）の結果に関しては、ベンチマークに用いたデータが少ないので、信頼性に欠けるかもとのこと。
🧸 実習：ホルムアルデヒドの励起状態の計算 ホルムアルデヒド（H2 CO}）について、
WebMO
を使って、下記で指定される計算手法で励起状態 (Excited State) 計算を実行してください。その際、必ず、振動解析の前に「構造最適化（Geometry Optimization） 」をおこなってください。Theory: RHF Basis Set: 6-31G(d) この分子の吸収スペクトルで観測される電子遷移のうち、最も光吸収の程度が大きいと予測される遷移 に対応する「吸収波長の値（nm 単位） 」と「振動子強度の値 」を調べてください。 上記が終わったら、計算手法（Theory / Basis Set）を変えて、計算結果、および、計算に掛かった時間を比較してみてください。
🍩 余談：円錐交差とは？ 分子の立体構造が変化しながら励起状態と基底状態の間のエネルギー差が段々と小さくなり、最終的に、そのエネルギー差がゼロになってしまう場合があります。 この「電子状態間のエネルギー差がゼロとなって互いに交差する」ことを円錐交差といいます。 Zgrablic, et al. (2011)
円錐交差する地点に辿り着いたとき、分子は、光を放出しないで基底状態に戻ります。 励起状態となった分子が「光の放出をともなわずに基底状態に戻る」ことを無輻射遷移 or 非断熱遷移といいます。非断熱遷移は、さまざまな機能性分子が多彩で巧みな機能を発現するための源です！ 分子の光化学反応を解析するときに、円錐交差の地点を探索し、その立体構造を調べたり、非断熱遷移の起こりやすさを調べたいと思うことがしばしばあります。
🍩 余談：エチレンの円錐交差点における分岐ベクトル たとえば、エチレン（CH2 =CH2 ）は、二重結合の回転にともなって S0 状態とS1 状態が近接し、円錐交差点に至ります。 このとき…分岐ベクトル $\mathbf{x}_1$：二重結合の回転方向 分岐ベクトル $\mathbf{x}_2$：ジラジカル構造 ↔ イオン構造の変化 Huix-Rotllant, et al. (2015)
🍩 余談：TD-DFT はエチレンの円錐交差点を正しく記述できない エチレンの円錐交差点における分岐ベクトル $\mathbf{x}_1$, $\mathbf{x}_2$ に沿ったエネルギー変化を見ると…SSR：State-interaction State-averaged spin-Restricted ensemble KS SF：Spin-Flip TD-DFT TD：TD-DFT TD-DFT（断熱近似）の場合、円錐交差の構造を適切に記述できていません。 Huix-Rotllant, et al. (2015)
🍩 余談：円錐交差点を適切に記述するには DFTスピン反転（spin-flip）TD-DFT アンサンブル（ensemble）DFT WFT
🤠 クイズ：電子励起状態の解析 量子化学計算に基づく電子励起状態の解析に関する次の文章について、適切なものをすべて選んでください。
Check
![Solov'yov, et al. (2012)]](https://gyazo.com/48f5e3d868003070ac0fbf8dd02f1e06.png)
