230228
https://gyazo.com/37fa642f89796e5e91a8ceaade04627e
資料へのリンク
https://gyazo.com/3704fe9705a09cf832ce8eddce572708
概要
日程： 2023 年 2 月 28 日（火）
会場： オンライン（Zoom）
目次
はじめに
学びの手引き
この資料について
資料については、講座が始まるまでに and/or 講座が終わってから、じっくりと読みながら学んでいただくことができるように準備しています
講座の終了後もできる限り公開を続けたいと思っているのですが、事情により、公開先を変更することがあるかも知れません
この講座の進め方
時間が限られてしまいますので、「ポイントのみを大まかに」という方針で説明することで、「量子化学計算をはじめようとするときに参考にしていただけそうな情報」をコンパクトにお伝えしたいと思っています
今回の講座は、量子化学計算をはじめるときの「足場」になればと思っています
星印（★）の項目は応用的な内容です
「はじめて量子化学計算に取り組もうとするとき」には、スキップしておいても良いように思います
疑問に思われたこと、詳細が知りたいと思われたこと、その他のこと、どんなことでも、Zoom の「Q&A ツール」からコメントをお送りください
代表的なツールやアプリの「実演」をおこないながら、「量子化学計算を実際におこなうときのイメージを掴んでいただくこと」も大切にしたいと思っています
演習課題・力試しテスト
学びの内容を振り返っていただくために演習課題・力試しテストも準備していますが、今回は十分に時間をとって解説することは難しいように思っています
●●という内容に対する演習課題・力試しテストの答えが気になる・解説してほしい、というご要望がありましたら、Zoom の「Q&A ツール」からお知らせください
プレゼンテーションモード
資料を共有する時には、プレゼンテーションモードという機能を使っています
ページの右側にある を押して、Start presentationを選ぶと、プレゼンテーションモードになります
この講座の講師
山本典史 (Yamamoto, Norifumi)
https://res.cloudinary.com/dfhjad2jo/image/upload/v1672805236/yamnor_o9bsoi.png
Profile
専門はコンピュータ化学、コンピュータを使って分子を解析しています
化学の学びを身近にすることにも興味を持っています
連絡先
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください
Links
この講座の到達目標
量子化学計算に取り組むときに必要となる基礎的な知識と基本的な技術について、説明できる
量子化学の知識や技術に基づいて、おもに有機化合物の化学的性質を解析するための指針について、説明できる
代表的な量子化学計算ソフトウェアである Gaussian と GAMESS の特徴と基本的な計算手順について、説明できる
代表的な量子化学計算ソフトウェアのインターフェースとして活用できる Web アプリ WebMO の特徴と基本的な使い方について、説明できる
この講座のスケジュール
table:schedule
開始 終了 時間 内容
10:30 12:00 90分 講義①
12:00 12:45 45分 昼食
12:45 14:15 90分 講義②
14:15 14:30 15分 休憩
14:30 16:00 90分 講義③
16:00 16:30 30分 質疑
16:30 終了
講義①
講義②
講義③
質疑
この講座の参考書
書籍版（3,850 円）
はじめて量子化学計算に取り組む人向けとして、読みやすい入門書です
Kindle 版（3,284円）
書籍版（3,520円）と Kindle 版（3,520円）
Gaussian に準拠した計算化学の入門書です
一般の書店では販売されていません
Gaussian の国内代理店などで購入できます
書籍版（10,000円）
GAMESS を使おうとするときに必要となる入力ファイル（計算の方法などを指定するもの）の具体的な書き方が丁寧に解説されているので、GAMESSを実際に使うときの大きな助けになります
書籍版（3,850円）
量子化学の基礎を学ぶ
量子化学計算の基礎となる「量子論」から、最近の研究開発分野で用いられている「密度汎関数法」のことまで、幅広く紹介されているバランスの良い教科書だと思います
上巻：書籍版（5,500円）と Kindle 版（5,225円）
下巻：書籍版（5,720円）と Kindle 版（5,434円）
書籍版（4,180円）と Kindle 版（3,971円）
計算化学をもっと深く学ぶ
書籍版（19,800円）
2023/3/29 に発売されるそうです
書籍版（4,049円）
書籍版（11,000円）
量子化学計算の概要
量子化学計算とは？
物質の様々な性質を特徴付けている微視的な物性諸量を知りたいと思うときに、対象としている物質の分子レベルでのふるまいを明らかにすることで、問題を解決するための様々な手がかりを与えてくれる便利な道具
機能性材料 を設計しようとしたり、新しい薬 を開発しようとするときに、研究・開発の現場では様々な課題に直面することもあります
たとえば、量子化学計算を活用すると、物質を構成している分子の形、安定性、分子同士の相互作用の強さ、化学反応の起こりやすさなどについて、ある程度の定量的な精度で信頼できる予測を理論的に導きだすことができます
分子シミュレーション手法の比較
table:計算対象と計算手法
対象 手法 有料アプリ 無料アプリ
分子系 第一原理 Gaussian GAMESS
材料系 第一原理 VASP Quantum Espresso
分子力場 Matlantis LAMMPS
生体系 分子力場 Amber Gromacs
各アプリの公式 Web サイト
第一原理計算のデータを深層学習したニューラルネットワークポテンシャル（NNP）力場を用いる
AmebrTools は無償配布
量子化学計算を使ってできること
分子の立体構造を予測できます
安定構造（最適化構造）
遷移状態の構造
構造異性体の存在比（ボルツマン分布）
分子の振動状態を予測できます
赤外スペクトル（振動スペクトル）
分子の励起状態を予測できます
吸収・発光スペクトル
円偏向性
分子の熱力学的な性質を予測できます
エンタルピー
エントロピー
定圧モル比熱
生成熱
分子の電子状態に基づいて様々な解析ができます
フロンティア軌道（HOMO や LUMO）
電子密度
部分電荷分布
分子の分極や溶媒効果を予測できます
双極子モーメント
溶媒和エネルギー
静電ポテンシャル
分子の化学反応を予測できます
反応経路
活性化エネルギー
力試しテスト
次に示すような「新しい薬を開発するケース」について、量子化学計算を用いることが適切だと思われるものをすべて選んでください。
x 薬物候補化合物として有力だと思われる幾つかの分子を合成するための準備として、合成経路の素反応のひとつについて、対応する遷移状態のエネルギーを知りたい。 x 薬物候補化合物として有力だと思われる幾つかの分子について、分子内でどのように電荷が分布しているか知りたい。 _ 100万個の化合物データベースの中から、薬物候補化合物としてふさわしい分子を1000個選びたい。 x 薬物候補化合物として有力だと思われる幾つかの分子の赤外吸収スペクトルについて、それぞれのピークがどのような振動モードに由来するのかを調べたい。 _ 薬物候補化合物として有力だと思われる幾つかの分子について、どのような結晶構造になっているのかを予測したい。 _ 薬物候補化合物として有力だと思われる幾つかの分子について、創薬のターゲットとなるタンパク質との相互作用の強さ（結合自由エネルギー）を予測したい。 量子化学計算の代表的なソフトウェア
業界標準の Gaussian
Pople は、いろいろあって、Gaussian の使用を拒否されてしまった 名前は、分子軌道を表現するために用いる ガウス関数（Gaussian） に由来
ガウス関数：$ \exp(-\alpha r^2)
スレーター関数：$ \exp(-\gamma r)
GaussView というアプリと併せて使うと、高機能な量子化学計算をお手軽に実行できる
https://i.gyazo.com/8ce4153a5a97629bd004dd4d55ccebca.png
計算化学者だけではなく、実験化学者にも Gaussian ユーザーが多い
分からないこと・困ったことがあるとき、適切なアドバイスをしてくれるユーザーがすぐに見つかる
インターネットの検索で有益な情報をたくさん得ることができる
最新版は Gaussian 16
ある敷地A内にある大学が、同敷地内の計算機センターの計算機上で外部の研究者に対してGaussianを利用させたい場合、コマーシャルサイトライセンスに加えてメインテナンスプログラムへの参加（有償）が必要
エネルギー計算 や 構造最適化 の収束性が良い
yamnor.icon 山本は、他の量子化学計算プログラムを使うときにも、構造最適化だけは Gaussian で実行する ことが多いです
Gaussian の入力ファイル
例：エチレン（$ \mathrm{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
学術＆商用で 無償利用可
yamnor.icon 商用でも無償という点は大きいと思います
オープンソース ではない
プログラムの再配布やソースコードの再利用は禁止
GAMESS-US を利用した研究成果を公表する場合、開発者らによる下記の 原著論文を引用する G. M. J. Barca, et al., J. Chem. Phys., Vol. 152, 154102 (2020)
アップデートが頻繁なので、バージョン番号ではなく、更新日で管理 されている
最先端で開発が進められている 様々な新しい手法 を実行できる
ソースコード を無償で入手することができる
オリジナルの機能を自分で追加する ことができる
量子化学計算に取り組む研究室では、派生バージョン を作っているところも多いのでは
yamnor.icon 山本も大学院生の頃、GAMESS を元にして、第一原理分子動力学（ab initio MD）のプログラムを作りました
GAMESS-US の入力ファイルの例
エチレン（$ \mathrm{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
イギリス版 GAMESS
イギリス国内 での学術用途に対しては 無償 で配布
イギリス国外 や 商用 での用途に対しては 有償 での配布
yamnor.icon 山本はまだ使ったことはありません
Firefly
モスクワ大学（ロシア） の Alexander Granovsky らが開発
無償で配布
1997年頃から GAMESS-US のソースコードを元にして開発が始まった
Intel 社製 CPU などで効率的に動作するように最適化されている
GAMESS-US から枝分かれしたときには PC-GAMESS という名前
2009年、Firefly という名前に変更
密度汎関数法を用いて量子化学計算をおこなった場合
GAMESS-US よりも計算にかかる時間が短い
最新バージョンは 2016年9月 に発表された 8.2.0
最近は更新されていない
yamnor.icon 山本はまだ使ったことはありません
この講座では
今回は主に、最も利用者が多いと思われる 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：静電ポテンシャルや双極子モーメントなど
演習課題：MolCalc を用いたエタノールの量子化学計算
MolCalc を用いてエタノール（$ \mathrm{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
計算方法と基底関数
量子化学計算の計算精度や計算コストは、主に、計算方法と基底関数の組み合わせで決まる
詳細はあとで
計算条件の表し方
計算方法の種類 と 基底関数の種類 を スラッシュ（/） で挟んで表記することが多い
たとえば...
B3LYP/6-31G(d)
MP2/cc-pVDZ
あるレベルで 構造最適化 した後、より高精度なレベルで エネルギー計算 をおこなった場合、２つの計算レベルを ダブルスラッシュ（//） で挟んで表記することがある
たとえば...
MP2/cc-pVDZ//B3LYP/6-31G(d)
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秒までの計算が可能
タイミングが悪いと混んでいる
上記のリンク先からアクセス
yamnor.icon 受講者の方には、アカウント情報をメールでお知らせしています
10分までの計算が可能
この講座の終了後、２週間は試していただくことができます
分子を SMILES で指定するときには、メニューから...
Lookup → Import by → SMILES
実演：WebMO を用いたエチレンの量子化学計算
https://scrapbox.io/files/63f99217e5f976001c763dc8.mov
演習：WebMO を用いたエタノールの量子化学計算
WebMO （講座用）を用いて、エタノール（$ \mathrm{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
量子化学計算でよく使うエネルギーの単位
1 Hartree = 627.5095 kcal/mol = 27.2114 eV
Hartree のことを、原子単位（atomic unit; a.u.）と書くこともある
量子化学計算のコスト
量子化学計算プログラム
Gaussian
GAMESS
学術 / 商用：0 円
GUI アプリ
GaussView
学術：約 16 万円
商用：約 30 万円
学術：1,200 USD（約 16 万円）
商用：3,000 USD（約 41 万円）
計算用サーバー
約 17 万円 / 年
力試しテスト
量子化学計算の代表的なソフトウェアに関する次の文章について、適切なものをすべて選んでください。
_ 量子化学計算のプログラムは数多くあるが、デファクトスタンダード（事実上の業界標準）だと言えるのは、アイオワ州立大学の John A. Pople らが中心に開発している GAMESS である。 x GAMESS-US は、オープンソースではなく、プログラムの再配布やソースコードの再利用は禁止されているが、学術・商用のどちらの用途でも無償で利用することができる。 x Gaussian は、有料の量子化学計算プログラムであり、GaussView というアプリと併せて使うことで、手軽に実行できる。また、エネルギー計算や構造最適化の収束性が良いことが知られている。 _ GAMESS には、現在、GAMESS-US、GAMESS-UK、Firefly という３つの異なるバージョンがあるが、全く同じ機能が実装されており、計算効率もほぼ同等である。 x ある分子について、Gaussian と GAMESS という異なるプログラムを用いて量子化学計算を実行しても、計算条件が等しければ、全エネルギーの値などについて、ほぼ等しい結果が得られる。 量子化学計算の代表的な計算方法
代表的な計算方法
電子相関を考慮しない計算方法
★ Møller-Plesset 摂動法
★ 結合クラスター法
★ 配置間相互作用（CI：Configuration Interaction）法
★ 多配置SCF法
★ 巨大分子系を扱う計算方法
★ フラグメント分子軌道（FMO：Fragment Molecular Orbital）法など
Hartree-Fock 法
ハートリー・フォック（Hartree-Fock; HF）法とは？
シュレディンガー方程式を解こうとするときに、電子同士の相互作用（電子相関）を厳密に取り扱うことができない
➡ なんとかして分子軌道を求めるために、何らかの「近似」（計算手法）が必要
HF法は、量子化学計算の出発点となる最も基本的な「近似」（計算手法）のひとつ
シュレディンガー方程式は３つ以上の多電子原子、あるいは分子を解析的に解くことができません（多体問題）。
yamnor.icon 追記
たとえば、原子と電子が１個ずつの「水素原子」はオーケーだけど、ヘリウム分子になると、原子１個と電子２個の計３体の系になるので、シュレディンガー方程式が解析的には解けなくなってしまいます
そこで使用するのが、ハートリー・フォック法（Hartree-Fock 法）です。
HF 法は、他の電子との相互作用を「平均化させたものとの相互作用」として近似することでこの問題を解決します
これを「一電子近似」と呼びます
https://gyazo.com/c67689935f786a7af9b12cc315beaa0c
まず、一電子近似で分子軌道を求めます
そこから、エネルギーが下がるように徐々に分子軌道を変化させ、もうこれ以上エネルギーが下がらないというところが実際の値に等しいと考えて、計算結果を出力させます
yamnor.icon 追記
「エネルギーが下がるように徐々に分子軌道を変化させる」ことを何度か繰り返すと、あるときに、「ある電子」と「他の電子との相互作用を平均化したもの（平均場）」との「つじつま」があうようになります
この「つじつま」があった状態のことを「自己無撞着な場（self-consistent field: SCF）」と呼びます
また、SCFを求めるために何度も繰り返し分子軌道を変化させる手続きのことを「SCF計算」と呼びます
Hartree-Fock 法の特徴
分子の立体構造や化学的性質について、多くの場合、十分な精度で予測できる
分子が最適構造にあるとき、その分子の全エネルギーの値を99%程度の精度で見積もることができる
しかし、化学反応を議論するためには、99%程度という精度でも不十分であり、より高精度な予測が必要となる
化学反応に伴う活性化エネルギー $ E_aを議論しようとする場合、$ E_aは 20 kJ/mol 程度（5 kcal/mol 程度）になることが多いので、この程度以下の誤差となることが望ましい
電子相関を考慮するように HF 法を拡張し、HF 法以上の計算精度を持つ手法のことを「post Hartree-Fock 法」と呼ぶ
post HF 法では、HF 法では無視してしまった「電子相関」の効果について、様々な工夫で考慮する
★ Hartree-Fock 方程式
シュレディンガー方程式
$ \left[ \frac{1}{2} \sum_i \nabla_i^2 - \sum_i \sum_\alpha \frac{Z_\alpha}{r_{i\alpha}} + \sum_{i} \sum_{j > i} \frac{1}{r_{ij}} \right] \Psi = E \Psi
３項目の $ \frac{1}{r_{ij}}は「電子-電子間の反発エネルギー」を表す
「ある電子は、他の電子の位置に依存しながら運動している」という描像
https://gyazo.com/c67689935f786a7af9b12cc315beaa0c
HF 方程式
$ \left[ \frac{1}{2} \sum_i \nabla_i^2 - \sum_i \sum_\alpha \frac{Z_\alpha}{r_{i\alpha}} + \sum_{i} V_i^{\mathrm HF} \right] \Psi = E \Psi
３項目の$ V_i^{\mathrm HF}は「ある電子-他の電子間の「平均」反発エネルギー」を表す
「ある電子は、他の電子の作る平均化された場の中で（他の電子と）独立に運動している」という描像
動的電子相関を考慮する計算方法
動的電子相関とは？
Hartree-Fock 法は、電子同士の相互作用を平均化して扱うという近似（平均場近似）をしている
このような電子同士の相互作用は、あとで登場する静的電子相関と区別するために、動的電子相関と呼ばれる
電子相関の大きさは「電子相関エネルギー」（$ E_\mathrm{corr}）と呼ばれ、一般に、シュレディンガー方程式の厳密解（$ E_\mathrm{exact}）と HF 計算で得られる値（$ E_\mathrm{HF}）との差分として定義される
$ 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) レベル
ただし、計算手法が高度になるほど、計算するためのコスト（計算に必要な時間やコンピュータの性能）も大きくなってしまうことがほとんどなので、目的に見合った精度とコストのバランスを考えることも大事になる
静的電子相関を考慮する計算方法
静的電子相関とは？
化学反応が起こる様子を追跡しようとすると、化学結合の解離や生成が起こるような領域では、複数の電子配置がエネルギー的に近接する場合があって、量子力学的共鳴に由来する電子間相互作用が顕著となる
このような相互作用のことを静的電子相関と呼ぶ
静的電子相関の具体例
典型例な静的電子相関の効果は、たとえば、水素分子のラジカル開裂にも現れる
https://gyazo.com/0b081d65040dfc7e3ab7680b49d5ec66
水素分子について、水素原子同士の結合を切断するために必要となるエネルギー（結合エネルギー）は、正確に解析してみると 4.75 eV になる
しかし、Hartree-Fock 法を用いて結合エネルギーを解析してみると、その値は 3.63 eV となり、正確な値とは大きく異なっている
このように定量性を欠いた結果となってしまうのは、水素分子の解離に伴って、この分子の基底電子状態と励起電子状態がエネルギー的に近くなり、２つの状態間の静的電子相関が顕著になるため
★ 静的電子相関を考慮する計算方法
HF 法は、電子状態を１つの電子配置のみで記述する近似法である
静的電子相関について、そのふるまいを正しく考慮するためには、電子状態を複数の電子配置で記述する多配置SCF（MCSCF：Multi Configurational SCF）法と呼ばれる方法が必要となる
Gaussian や GAMESS では、最も一般的な MCSCF 法である完全活性空間 SCF（CASSCF：Complete Active Space SCF）法をはじめとして、静的電子相関の寄与を考慮するためのさまざまな計算方法を使うことができる
★ 静的電子相関＋動的電子相関を扱う方法
多参照 SCF（MRSCF：Multi Reference SCF）法は、静的電子相関を記述する MCSCF 法で求めた波動関数に対して動的電子相関を考慮する方法である
GAMESS では、摂動論に基づくMRMP（Multi Reference Møller-Plesset）法などを利用することで、金属錯体の励起状態など、静的・動的電子相関が絡む複雑な電子状態に関しても定量的精度で各種物性を求めることができる
密度汎関数法
密度汎関数法（DFT：Density Functional Theory）とは？
近年、研究開発の最前線で大活躍している計算手法
電子密度を試行関数とする多電子状態問題の近似解法
（動的）電子相関の寄与について、 Hartree-Fock 法とほぼ同程度の計算コストで考慮することができる
Hartree-Fock 法 などは波動関数を試行関数とする方法であるが、DFT は電子密度を試行関数とする方法であり、その基盤となる考え方が大きく異なっている
HF 法の基礎は HF 方程式だが、DFT の基礎はコーン・シャム（Kohn-Sham; KS）方程式
密度汎関数法が得意なこと
DFT は、電子間相互作用を含む汎関数を導入することで、HF 方程式とほぼ同様な KS 方程式を解くだけで、動的電子相関の寄与を考慮することが可能
動的電子相関の寄与は考慮できるが、基本的には、静的電子相関の寄与は考慮できない
したがって、利用する汎関数が適切であれば、HF 法を使うくらいの安価な計算コストで、MP2 摂動法レベルの電子相関効果を含む精度の解析結果を得ることができる
このような魅力から DFT は、HF 法が苦手とする金属含有系をはじめとして、幅広い物質群の応用計算に利用されている
密度汎関数法が苦手なこと
一方で DFT の弱点として、長距離電子相関の見積もりは不得意としている
分子同士に働く相互作用のうち、ファンデルワールス（van der Waals）相互作用（分散力）を過小評価してしまい、分子集合体系の安定性などを正しく評価できない
長距離領域での電子間相互作用を補正する「長距離補正（Long-range Correction; LC）法」などの改良が続けられている
★ KS 方程式とは？
HF 方程式
$ \left[ \frac{1}{2} \nabla^2 - \sum_\alpha \frac{Z_\alpha}{r_{\alpha}} + V^\mathrm{HF} \right] \psi_k = \epsilon_k \psi_k
３項目の$ V^\mathrm{HF}は「ある電子-他の電子間の「平均」反発エネルギー」を表す
「ある電子は、他の電子の作る平均化された場の中で（他の電子と）独立に運動している」という描像
KS 方程式
$ \left[ \frac{1}{2} \nabla^2 - \sum_\alpha \frac{Z_\alpha}{r_{\alpha}} + V^\mathrm{eff} \right] \psi_k = \epsilon_k \psi_k
３項目の $ V^\mathrm{eff}は、「外場有効ポテンシャル」
$ V^\mathrm{eff} = V^\mathrm{Coulomb} + V^\mathrm{XC}
$ V^\mathrm{Coulomb}：電子-電子間の古典的な静電（クーロン）相互作用
$ V^\mathrm{XC}：非古典的な相互作用であり、「（交換相関）汎関数」と呼ばれる
この「（交換相関）汎関数」は、DFT を特徴付ける重要なもの
KS 方程式とHF 方程式の違いは、３項目が「平均場$ V^\mathrm{HF}」か「外場有効ポテンシャル$ V^\mathrm{eff}」かのみ
つまり、KS 方程式は、HF 方程式と同程度の時間で計算できる
HF 方程式は「シュレディンガー方程式の近似式」であるのに対して、KS 方程式は「シュレディンガー方程式と同等」である
つまり、もし KS 方程式が解けるのならば、計算時間は HF 法と同程度であるにも関わらず、厳密な解が得られるはず
HF 方程式中の「平均場$ V^\mathrm{HF}」は「既知の関数」であるが、KS 方程式中の「交換相関汎関数$ V^\mathrm{XC}」は厳密な形が不明な「未知の関数」である
$ V^\mathrm{XC}の厳密な形が分からないので、人為的に仮定された汎関数を使う
この「人為的に仮定された」汎関数の質が、DFT の精度を決定することになる
密度汎関数法で用いる「汎関数」
汎関数とは？
DFT の計算精度と信頼性は、電子密度と系のエネルギーを関係付ける汎関数の種類に依存する
適切な汎関数を選べば、良い計算精度が得られる
汎関数の種類
基本的な汎関数
BLYP （Becke-Lee-Yang-Parr）
PBE（Perdew-Burke-Ernzerhof）
Hartree-Fock 交換項と組み合わせたハイブリッド型の汎関数
B3LYP
以前は実質上の業界標準（デファクトスタンダード）であるかのように多用されていました
yamnor.icon 山本は、今でも、汎関数に迷ったときの第一候補かなと思っています
長距離電子相関の見積もりが不得意であることが知られるようになった
長距離補正を含む「CAM-B3LYP」、分散力補正を含む「B3LYP-D3」などが開発されている
PBE0
PBEPBE
長距離補正を含む汎関数
CAM-B3LYP
ωB97シリーズ（ωB97, ωB97-X, ωB97-XD）
LC-BOP
分散力補正を含む汎関数
B3LYP-D3
ωB97-XD
yamnor.icon 山本は分子の励起状態を計算するときに、これをよく使います
ミネソタシリーズ（M06-2X, X11 など）
力試しテスト
Hartree-Fock (HF) 法に関する次の文章について、適切なものをすべて選んでください。
x HF 法は、「ある電子について、他の電子の作る平均化された場の（平均場）中で、他の電子とは独立に運動している」という描像で近似している。 _ HF 法は、量子力学的共鳴に由来する「静的電子相関」について、定量的な精度で十分に考慮することができる計算方法である。 _ HF 法は、電子同士の相互作用に由来する「動的電子相関」について、定量的な精度で十分に考慮することができる計算方法である。 x HF 法は、分子の立体構造や化学的性質について、多くの場合、十分な精度で予測できる。たとえば、全エネルギーの値についても、99%程度の精度で見積もることができる。 _ HF 法は、化学反応に伴う活性化エネルギーについて定量的な精度で議論したい場合でも、十分に信頼できる結果が得られる計算方法である。 力試しテスト
密度汎関数理論（DFT）法に関する次の文章について、適切なものをすべて選んでください。
x DFT は、電子同士の相互作用に由来する「動的電子相関」について、定量的な精度で十分に考慮することができる計算方法である。 _ DFT は、量子力学的共鳴に由来する「静的電子相関」について、定量的な精度で十分に考慮することができる計算方法である。 _ DFT の汎関数として代表的な B3LYP を使った場合、分子集合体系の安定性などについても、定量的な精度で正しく評価できる。 _ DFT の計算精度と信頼性は、どのような汎関数を用いた場合でも、ほぼ同様である。 _ DFT は、HF と同様に、波動関数を試行関数とする方法であり、その基礎となる Kohn-Sham 方程式は、シュレディンガー方程式を近似したものであると言える。 量子化学計算の代表的な基底関数
基底関数とは？
任意の波動関数 $ \Psi を有限個の既知の関数 $ \phi_1, \phi_2, ..., \phi_n の線形結合で表すという近似（線形近似）を考える
$ \Psi = c_1 \phi_1 + c_2 \phi_2 + ... + c_n \phi_n = \sum_i^n c_i \phi_i
この既知の関数 $ \phi_1, \phi_2, ..., \phi_n のことを 基底関数 と呼ぶ
展開係数 $ c_1, c_2, ..., c_nの値を変えると、波動関数 $ \Psiが変わる
量子化学計算では、あらかじめ適切にパラメータを設定した基底関数 $ \{\phi_n\}を使って、展開係数 $ \{c_n\}の値について、最適なものを探そうとする
つまり、計算の途中で...
基底関数 $ \{\phi_n\}は変えずに、
展開係数 $ \{c_n\}の値だけを変える。
https://i.gyazo.com/a4e3310a0265db15ba8d23dbb445192c.png
基底関数には様々な種類があり、どのような分子を解析したいのか、何を調べたいと思っているのか、どのような計算方法を用いるのか、それぞれに応じて、適切に選択する必要がある
基底関数の個数が多いほど、電子分布を柔軟に表現できるので、計算精度が高くなる
計算精度が高くなればなるほど、計算コストも高くなってしまう
計算精度と計算コストのバランスを考えて、適切に選ぶ
基底関数の「数」と分子軌道の「数」
$ 2N個の電子を持つ分子について、$ L個の基底関数を用いて量子化学計算をすると、通常...
占有軌道（電子が詰まった軌道）：$ N個
$ N個の占有軌道のうち、最もエネルギーが高いものが HOMO（Highest Occupied Molecular Orbital）
非占有軌道（電子が詰まっていない軌道）：$ L-N個
$ L-N個の非占有軌道のうち、最もエネルギーが低いものが LUMO（Lowest Unoccupied Molecular Orbital）
基底関数の個数を増やしていくと、非占有軌道の数が多くなる
計算精度は高くなるが、「化学的な意味を与えにくいもの」が多く含まれるようになる
https://i.gyazo.com/122f5015a31b0d670b4fc9bbda78f820.png
水素原子の場合
水素分子を表す波動関数 $ \Psi は、たとえば、２個の水素原子 $ \mathrm{H}_aと$ \mathrm{H}_bの中心に置いた基底関数$ \phi_aと$ \phi_bの線形結合として、次のように表すことができる
$ \Psi = c_a \, \phi_a + c_b \, \phi_b
https://gyazo.com/0857066e42c0815178e589a5802444d1
https://gyazo.com/dbb67c916e73c77f44cf84a176b59e55
https://gyazo.com/d3e4838821594b53661bac97ddfcb60d
水素原子については、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は、関数の「広がりの程度」を表す変数
水素原子の厳密な解析解（= 電子の分布を適切に表すことができる）なので、計算精度は高い
ただし、数値計算の計算効率は悪い
yamnor.icon 図にあるように、原子の中心（$ r = 0）で値が不連続になるから、数値的に積分することが難しい
計算効率が悪いので、ふつうは、量子化学計算では使わない
量子化学計算では、スレーター型関数の代わりに、$ \exp(-\alpha r^2)型の関数を基底関数として使う
$ \phi_a = \exp(-\alpha r^2)
このような関数を ガウス型関数 と呼ぶ
https://gyazo.com/66a23423bf1cc335cf3e817cdf51f374
係数 $ \alphaは、関数の「広がりの程度」を表す変数
量子化学計算の前に、あらかじめ設定されている定数（計算の途中で変化させない）
ガウス型関数は、原子軌道に対する近似解なので、計算精度は高くはない
ただし、数値計算の計算効率が良い
yamnor.icon 上の図にあるように、原子の中心（$ r = 0）でも値が連続になるから、数値積分が簡単
計算精度の悪さは、拡がりの程度を変えたガウス型の基底関数をいくつか線形結合させることでカバーする
いくつかのガウス型関数を線形結合させたものを短縮型ガウス基底とよぶ
たとえば...
$ \phi_\mathrm{赤} = c_\mathrm{緑} \, \chi_\mathrm{緑} + c_\mathrm{紫}\, \chi_\mathrm{紫} + c_\mathrm{黄}\, \chi_\mathrm{黄}
$ = 0.5 e^{-0.5 r^2} + 0.3 e^{-1.0 r^2} + 0.2 e^{-0.3 r^2}
https://gyazo.com/40534f7c97701c1c01c4e82ecef357c2
赤い実線で示す短縮型ガウス関数は、青い点線で示すスレーター型関数のかたちを、まあまあ良く再現している
基底関数の種類①
最小基底系（STO-nG基底系）
例
STO-3G
３種類のガウス型関数（3G）を足し合わせて、スレーター型軌道（STO）に似せたもの
STO-4G
STO-5G
計算精度の比較
STO-3G < STO-4G < STO-5G
計算コストは低いが、計算精度も低いので、実験結果と比較するなどの実用的な解析に用いることはあまりない
★ 実際の配置
例
第１周期原子
１つの s 軌道のみで構成される
表記：[1s]
第２周期原子
2つの s 軌道と 1つの p 軌道で構成される
表記：[2s1p]
基底関数の種類②
Pople 系
例
3-21G
6-31G
6-311G
計算精度の比較
3-21G < 6-31G < 6-311G
最小基底系に比べると、より柔軟に電子のふるまいを記述できるため、計算精度が改善される
原子価軌道に対して、拡がりの異なる複数の関数を組み合わせて表現する
一般的な有機化合物に対しては、基盤となる 6-31G に 分極関数 を加えた 6-31G(d) や 6-31G(d,p) などが、計算効率が良い・実用的な精度をもつ基底関数としてよく用いられている
★ 実際の配置
6-31Gの場合
１つの占有軌道に対して...
内殻軌道：６個のガウス型を足し合わせた（短縮した）１種類の関数で記述
原子価軌道：２個のガウス型を短縮したもの＋別の１個のガウス型を組み合わせて、計２種類の関数で記述
２種類の関数で原子価軌道を表すものを double zeta 型とよぶ
第１周期原子
原子価軌道：[2s]
たとえば、水素（H）の場合
1s 軌道 → 2 個
基底関数は合計で 2 個
第２周期原子
内殻軌道：[1s]
原子価軌道：[2s2p]
たとえば、炭素原子（C）の場合
1s軌道 → 1 個
2s 軌道 → 2 個
2p 軌道 → 2 × 3 = 6個
基底関数は合計で 9 個
基底関数の種類③
高精度計算に適した基底関数（Correlation-Consistent 基底系）
計算精度の比較
cc-pVDZ < cc-pVTZ < cc-pVQZ
応用研究の中では、double zeta 精度の cc-pVDZ が頻繁に用いられる
double zeta 精度では、電子相関効果の取り込みが不十分である場合もある
化学反応などにおいて 10 kcal/mol 程度のエネルギー差を定量的に議論したいときには、cc-pVTZ 程度（triple zeta 精度）の基底関数を CCSD(T) レベルの計算手法と組み合わせて用いる
分極関数とは？
分子の複雑な構造を適切に記述したい場合には、各原子上に配置した基底関数に 分極関数（polarization function） を追加することで、電子分布が柔軟に変形する自由度を与えることができる
たとえば...
p 軌道をもつ水素以外の原子に対しては、d 型の関数を追加する
http://4.bp.blogspot.com/-9odFiZctEe0/VX1_xMr2R6I/AAAAAAAADi4/uLpgZkL2MU0/s1600/Fig2.7.8.jpg
球対称な s 軌道のみを持つ水素原子に対しては、p 型の関数を追加する
http://2.bp.blogspot.com/-nQJzi068Llo/VX1_kBl-cZI/AAAAAAAADiw/aE4Ln8gsYFM/s320/Fig2.7.7.jpg
表記方法 / 指定方法
Pople 系に対して、分子内にある水素以外の原子に分極関数を追加する場合
3-21G(d) or 3-21G*
6-31G(d) or 6-31G*
Pople 系に対して、水素以外の原子に分極関数を追加した上で、水素原子にも分極関数を追加する場合
3-21G(d,p) or 3-21G**
6-31G(d,p) or 6-31G**
cc-pVnZ 基底関数（$ n = D, T, Q, 5, 6）は、水素原子および水素以外の原子のどちらにも、全ての原子に分極関数が追加される
分散関数とは？
負電荷を帯びた分子は、中性分子と比較して、電子分布が広がっている場合がある
このような場合、物性諸量を正しく見積もるためには、電子分布の広がりを記述するために 分散関数（diffuse function） を加えることができる
表記方法 / 指定方法
Pople 系に対して、分子内にある水素以外の原子に分散関数を追加する]場合
3-21+G(d) or 3-21+G*
6-31+G(d) or 6-31+G*
Pople 系に対して、水素以外の原子に分散関数を追加した上で、水素原子にも分散関数を追加する場合
3-21++G(d,p)
6-31++G(d,p)
cc-pVnZ 基底関数（$ n = D, T, Q, 5, 6）に対して、水素原子および水素以外の原子のどちらにも、全ての原子に分散関数を追加する場合
aug-cc-pVDZ
aug-cc-pVTZ
オススメの計算手法と基底関数は？
yamnor.icon は、一般的な有機化合物の基底状態の物性を調べようとする場合、まずは B3LYP/6-31G(d) や B3LYP/6-31G(d,p) を使うことが多いです
原子間の結合距離を定量的に議論するには、最低でも double zeta 精度（6-31G）が必要
原子間の結合角度を定量的に議論するには、分極関数（6-31G(d)）が必要
★ 周期性を持つ分子系は？
固体表面や結晶などの周期性を持つ分子系に興味がある場合には、原子を展開中心とするガウス型基底を用いてその電子状態を解析しようとすると難しい
したがって、一般的には、量子化学計算が得意な解析可能な対象は、主として孤立分子系に限られる
周期性を持つ分子系に対しては、平面波を基底関数とする計算化学プログラム（たとえば VASP など）の方が助けになる
力試しテスト
基底関数に関する次の文章について、適切なものをすべて選んでください。
_ 量子化学計算では、ガウス型関数の代わりに、スレーター型関数を基底関数として使う。なぜならば、スレーター型関数は、計算精度が高いうえに、計算効率も良いからである。 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) を用いた方がよい。 量子化学計算に基づく構造最適化
分子の安定構造の探索
各原子核に働く力（エネルギー勾配 $ \frac{\partial E}{\partial Q}）の方向にカタチをちょっとずつ変化（$ x_\mathrm{new} = x_\mathrm{old} - \alpha \frac{\partial E}{\partial Q}）させることで、分子はより安定となる（ = エネルギーが小さくなる）
$ E(x_\mathrm{new}) < E(x_\mathrm{old})
https://gyazo.com/6fad80947159d295366c18cfb474210c
「カタチをちょっとずつ変化させる」を繰り返すことで、分子の立体構造を半自動的に最適化できる
実際には、多くの量子化学計算プログラムは、エネルギー勾配の方向に動かすという単純な方法（勾配降下法）ではなく、もっと効率的なアルゴリズムを使っている
Gaussian では、デフォルトでは、GEDIISというアルゴリズムを用いる
実演：Gaussian を用いた水分子の安定構造の探索
https://scrapbox.io/files/63f9a0bcd06e52001c9d4a71.mov
最適化される様子を確かめるため、HOH 結合角を 150° 程度にしたものを初期構造とした
分子の遷移状態の探索
ある化学反応について、反応速度を議論するときには、その反応に伴う遷移状態を調べる必要がある
たとえば、ある反応の速度定数 $ kについて、始状態と遷移状態とのエネルギー差、つまり、活性化エネルギー $ \Delta G^{\ddagger}が分かっていれば、Eyring-Polanyi 式から...
$ k = \chi \frac{k_\mathrm{B}T}{h}\exp\left(-\frac{\Delta G^{\ddagger}}{RT}\right)
$ \chi：遷移状態において、反応座標に沿った分子振動が一往復したときに反応する確率（透過係数）
例：Diels-Alder 反応
https://gyazo.com/6245b572e337bc0273b0b260af89a429
https://i.gyazo.com/a53e519a0a2b7c54971c9b3bd1ef1871.png
多くの量子化学計算プログラムには、遷移状態を探索する機能が備わっている
単純に「構造最適化の逆」（$ x_\mathrm{new} = x_\mathrm{old} + \alpha \frac{\partial E}{\partial Q}）をやっても、うまくいかない
https://gyazo.com/6fad80947159d295366c18cfb474210c
エネルギーの二次微分（$ \frac{\partial^2 E}{\partial Q^2}）を要素にもつ Hessian 行列 の固有値（力の定数 or 曲率：$ \chi）を求めて...
Hessian 行列の負の固有値を持つ固有ベクトルの方向に少しずつ動かすという方法がある
→ Eigenvector Following 法
調和近似（$ E = \kappa Q^2）からの歪みが大きな方向に少しずつ動かすという方法がある
Hessian 行列を計算するためのコストがまあまあ大きい
構造最適化よりも時間がかかる
実演：Gaussian を用いた水分子の遷移状態の探索
https://scrapbox.io/files/63f9c6f1a1e4d5001b464477.mov
遷移状態に近い（と予想される）構造から出発する
今回の場合、H-O-H 結合角を大きくする
yamnor.icon 見つけようとしている遷移状態は直線型だと考えられるから
遷移状態探索の設定方法
Job Type
https://i.gyazo.com/c1ba79cb670331d2f93de606d041500f.png
Optimizationを選択
Optimize to a TS (Berry)
Force Constants Calculate at First Point
遷移状態で振動解析すると...
値が負（虚数）となる振動数がある
https://i.gyazo.com/862a1da82d62e9ec3de7e7916cd77573.png
安定構造では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
局所安定状態（下に凸の放物線）であり、エネルギーの座標に関する二次微分の値が常に正である
$ \frac{\partial^2 E}{\partial Q^2} > 0
振動解析すると、すべての振動モードの振動数が正の値になる
https://gyazo.com/e949b48ac2e562a22b8de1963cf0c3ec
遷移状態では...
エネルギーの座標に関する一次微分（エネルギー勾配）が常にゼロである
$ \frac{\partial E}{\partial Q} = 0
特定の変位方向に対して局所不安定状態（上に凸の放物線）であり、その変位方向にはエネルギーの座標に関する二次微分の値が負となる
$ \frac{\partial^2 E}{\partial Q^2} < 0
振動解析すると、１つの振動モードの振動数のみが負の値になる
https://gyazo.com/1f0cfca95edab362017ddd028edbef68
演習課題：水分子の構造最適化
WebMO（講座用）を用いて、水（$ \mathrm{H_2O}）について、下記の計算方法と基底関数の中から好きな組み合わせを選んで構造最適化（Geometry Optpmization）し、実験値と比較してみてください。 計算方法
Hartree-Fock (HF)
DFT / B3LYP
Moller Plesset 2 (MP2)
基底関数
3-21G
6-31G(d)
6-311+G(d,p)
実験値
O-H 間の結合長
0.9575 Å
H-O-H 間の結合角
104.51°
計算方法と基底関数の組み合わせ（水分子の構造最適化の例）
table:水分子の構造最適化
計算方法 基底関数 結合長 誤差 結合角 誤差
（Å） （°）
HF 3-21G 0.967 +0.01 107.69 +3.18
6-31G(d) 0.947 -0.01 105.51 +1.00
6-311+G(d,p) 0.941 -0.02 106.22 +1.71
B3LYP 3-21G 0.997 +0.04 103.91 -0.60
6-31G(d) 0.969 +0.01 103.64 -0.87
6-311+G(d,p) 0.962 +0.01 105.05 +0.54
MP2 3-21G 0.989 +0.03 105.16 +0.65
6-31G(d) 0.969 +0.01 103.97 -0.54
6-311+G(d,p) 0.960 +0.00 103.37 -1.14
【演習課題３】構造最適化
構造最適化に関する次の文章について、適切なものを全て選んでください。
_ 最適構造では、エネルギーの座標に関する一次微分（$ \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 行列に基づく方法を使うので、安定構造の探索よりも計算コストがかかる。 量子化学計算に基づく振動解析
振動解析とは？
分子固有の振動数と赤外強度を理論的に計算する方法
実験的に測定された赤外スペクトルの解釈を補助する手段として活用されている
各振動モードに対するラマン強度を計算することも可能
赤外吸収スペクトルとは？
水分子の赤外吸収スペクトル（実験結果）
https://gyazo.com/246c01937cec7d08b6c54b43d41e6f3e
分子振動とは？
水（$ \mathrm{H_2O}）の変角振動
https://gyazo.com/e4d9bb89f0f61f0533b9532d0e024287
水（$ \mathrm{H_2O}）の対称伸縮振動
https://gyazo.com/145fab5f0795498d91045fd285f6b556
水（$ \mathrm{H_2O}）の逆対称伸縮振動
https://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
★ 分子振動を表す式（古典力学の場合）
原子１と２から成る２原子分子を考える
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
演習課題：水分子の振動解析
WebMO（講座用）を用いて、水分子（$ \mathrm{H_2O}）について、好きな「計算方法 & 基底関数 (Optimize + Vib Freq)」の異なる組み合わせを用いて「構造最適化」と「振動数解析」 をおこない、得られた計算結果を実験値と比較してみてください。 実験結果
変角振動：1595 $ \mathrm{cm}^{-1}
対称伸縮振動：3657 $ \mathrm{cm}^{-1}
逆対称伸縮振動：3756 $ \mathrm{cm}^{-1}
計算方法と基底関数の組み合わせ（水分子の振動解析の例）
table:水分子の振動解析
計算方法 基底関数 変角 誤差 対称伸縮 誤差 逆対称伸縮 誤差
(cm-1) (cm-1) (cm-1)
HF 3-21G 1799 204 3812 155 3946 190
6-31G(d) 1826 231 4071 414 4189 433
6-311+G(d,p) 1726 131 4141 484 4244 488
B3LYP 3-21G 1693 98 3411 -246 3551 -205
6-31G(d) 1713 118 3724 67 3846 90
6-311+G(d,p) 1602 7 3813 156 3919 163
MP2 3-21G 1722 127 3504 -153 3663 -93
6-31G(d) 1735 140 3775 118 3917 161
6-311+G(d,p) 1630 35 3880 223 3999 243
分子振動の非調和性
高い精度（MP2/cc-pVTZ）で計算した場合
https://gyazo.com/e1d212dc0d548ac021fe5ffb7701964a
MP2摂動法を用いることで、実験値にかなり近くなる
高波数の振動モード（対称伸縮＆逆対称伸縮）では、誤差が大きい
分子振動を調和振動子だと見なす近似（調和近似）による誤差（非調和性）の影響だと考えられる https://gyazo.com/dcff54c2ebf09e24e2db0b7bfa25ff73
非調和性の問題を手軽に解決する場合には、振動数に適当な数字（スケーリング因子）を掛けて補正する
事例：
HF/6-31G(d) の場合 → 0.903
B3LYP/6-31G(d,p) の場合 → 0.961
https://i.gyazo.com/ec1b28045886c70d80aa644f4f88a593.png
★ 非調和振動数解析
分子振動をより定量的精度で解析したい場合には、調和近似を越えた取り扱いが必
GAMESS では、非調和性を考慮した振動解析として、Vibrational SCF (VSCF) 法を利用することができる
VSCF法で計算した場合
https://gyazo.com/d66cdc593c263cd07ac4f263e4ba3715
量子化学計算に基づく化学反応の解析
化学反応を理解する上で重要な考え方は？
化学平衡論
反応前後で変化する自由エネルギー（エンタルピー）の大きさを知ることによって、反応前後の状態に対する分布比が予測できる
反応速度論
遷移状態における自由エネルギー（エンタルピー）の高さを知ることによって、化学反応の速度定数が予測できる
化学平衡論に基づくと…
ある分子の集団が熱平衡にあるとき、1 mol あたりのなかで、ある特定の状態$ iを取り得る確率は...
$ P_i \propto \exp\left(-\frac{\epsilon_i}{RT} \right)
$ \epsilon_i：ある状態$ iにおけるエネルギー
$ R：気体定数（$ = k_\mathrm{B} N_\mathrm{A} = 1.987\,206 \; \mathrm{cal} \; \mathrm{K}^{-1} \; \mathrm{mol}^{-1}）
$ k_\mathrm{B}：ボルツマン定数（$ = 1.380649 \times 10^{−23} \; \mathrm{J} \; \mathrm{K}^{−1}）
$ N_\mathrm{A}：アボガドロ定数（$ = 6.02214076 \times 10^{23} \; \mathrm{mol}^{−1}）
同様に、1 mol あたりのなかで、ある分子が２つの状態$ iと$ jであるときの存在比は...
$ \frac{P_i}{P_j} = \exp\left(-\frac{\Delta \epsilon_{ij}}{RT}\right)
$ \Delta \epsilon_{ij}：エネルギー差（$ = \epsilon_i - \epsilon_j）
$ \epsilon_i：ある状態$ iにおけるエネルギー
$ \epsilon_j：ある状態$ jにおけるエネルギー
たとえば、ある分子について、２つの状態$ iと$ jのエネルギー差が$ \Delta \epsilon_{ij} = 5 \; \mathrm{kcal} \; \mathrm{mol}^{-1}である場合、温度$ T = 300 \; \mathrm{K}における存在比は...
$ \frac{P_i}{P_j} = \exp \left(-\frac{5.0 \times 10^3 \;\mathrm{cal\;mol^{-1}}}{1.987\,206 \; \mathrm{cal\; K^{-1} \; mol^{-1}}\times 300 \;\mathrm{K}} \right) = 2.278 \times 10^{-4}
この場合、状態$ iは$ jに対して、0.02 % 程度しか存在しなと予想される。
反応速度論（アレニウス則）に基づくと…
ある素反応過程について、始状態と遷移状態のエネルギー差（エンタルピー差）が$ \Delta E^{\ddagger}であるとき、反応速度定数は...
$ k = A \exp\left(-\frac{\Delta E^{\ddagger}}{RT}\right)
$ A：頻度因子（前指数因子）
たとえば、ある反応について、始状態と遷移状態のエネルギー差が$ \Delta E^{\ddagger} = 5 \; \mathrm{kcal} \; \mathrm{mol}^{-1}であり、頻度因子が$ A = 1.0 \times 10^{10} \; \mathrm{sec}^{-1}である場合、温度$ T = 300 \; \mathrm{K}における反応速度定数は...
$ k = 1.0 \times 10^{10} \; \mathrm{sec}^{-1} \; \times\exp \left(-\frac{5.0 \times 10^3 \; \mathrm{cal\;mol^{-1}}}{1.987\,206 \; \mathrm{cal\; K^{-1} \; mol^{-1}} \times 300 \; \mathrm{K}} \right) = 2.278 \times 10^{6} \; \mathrm{sec}^{-1}
反応速度論（遷移状態理論）に基づくと…
ある素反応過程について、始状態と遷移状態のエネルギー差（エンタルピー差）が$ \Delta E^{\ddagger}であるとき、反応速度定数は...
$ k = \frac{k_\mathrm{B} T}{hc} \exp\left(-\frac{\Delta E^{\ddagger}}{RT}\right)
$ k_\mathrm{B}：ボルツマン定数
$ h：プランク定数
$ c：濃度
実演：Gaussian を用いたブタジエンの cis/trans 異性化反応の解析
https://scrapbox.io/files/63faad6c14c7b1001bf5afac.mov
演習課題：ブタジエンの 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}となります。量子化学計算を実行する際には、基底関数は 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
問４
ブタジエンの trans 体と cis 体のエネルギー差 $ \Delta E_\mathrm{trans-cis} = \epsilon_\mathrm{trans} - \epsilon_\mathrm{cis}を$ \mathrm{kcal} \; \mathrm{mol}^{-1}単位で答えてください
$ \Delta E_\mathrm{trans-cis}= -3.547 $ \mathrm{kcal \; mol^{-1}}
問５
ブタジエンの trans 体と cis 体の存在比$ \frac{P_\mathrm{trans}}{P_\mathrm{cis}}の値を答えてください
$ \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 体と遷移状態のエネルギー差$ \Delta E_\mathrm{TS-trans} = \epsilon_\mathrm{TS} - \epsilon_\mathrm{trans}を$ \mathrm{kcal} \; \mathrm{mol}^{-1}単位で答えてください
$ \Delta E_\mathrm{trans-TS} = 7.561 $ \mathrm{kcal \; mol^{-1}}
問７
ブタジエンの cis 体と遷移状態のエネルギー差$ \Delta E_\mathrm{TS-cis} = \epsilon_\mathrm{TS} - \epsilon_\mathrm{cis}を$ \mathrm{kcal} \; \mathrm{mol}^{-1}単位で答えてください
$ \Delta E_\mathrm{cis-TS} = 4.014 $ \mathrm{kcal \; mol^{-1}}
問８
ブタジエンの trans → cis 異性化反応に対する反応速度定数$ k_{\mathrm{trans}\rightarrow\mathrm{cis}}を$ \mathrm{sec}^{-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}単位で答えてください
$ 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://gyazo.com/0b9a593af6e116be42aa973d6e203c8b
$ E_{励起} - E_{基底} = \Delta E = h c / \lambda_{abs}に相当する光のエネルギーを分子が吸収し、基底 → 励起状態への電子遷移が起こる
$ h = 4.1357 $ \times 10^{-15} (eV • s)
$ h：Planck 定数
$ c：光の速度
$ c = 2.9979 $ \times 10^8 m/s
$ h c = 1.2398 $ \times 10^{-6} eV • m = 1240 eV • nm
共役分子系が光を吸収すると?
$ \pi軌道にある電子が、$ \pi^{*}軌道へ遷移する
https://gyazo.com/c6e786be2c17fb1bd8b13f77756d8428
$ \pi\pi^{*}遷移と呼ぶ
$ \Delta E_{\pi\pi^{*}} = \frac{1240}{\lambda_\mathrm{abs}}
$ \Delta E_{\pi\pi^{*}}：遷移エネルギー［ev］
$ \lambda_\mathrm{abs}：光の吸収波長［nm］
鎖状共役系の吸収スペクトルは？
分子鎖が長くなると、吸収波長が長くなる
https://i.gyazo.com/a19826a07e1c8bc43341c359dab0f4e6.png
分子の吸光・発光と基底・励起状態
https://gyazo.com/cc8fd20bf491678066145f87d47dc9a6
吸光
基底（$ \mathrm{S_0}）状態の最適構造 → 励起（$ \mathrm{S_n}）状態
つまり、量子化学計算に取り組むときに、分子の立体構造は、基底状態で最適化する必要がある
発光
励起（$ \mathrm{S_n}）状態の最適構造 → 基底（$ \mathrm{S_0}）状態
つまり、量子化学計算に取り組むときに、分子の立体構造は、励起状態で最適化する必要がある
励起状態の計算方法
時間依存密度汎関数理論（Time-Dependent DFT）法。励起状態の計算方法として、最も普及している
一電子励起配置間相互作用法
Symmetry Adapted Cluster/Configuration Interaction (SAC-CI) 法
完全活性空間 (Complete Active Space) SCF 法
実演：ホルムアルデヒドの電子励起状態の解析
https://scrapbox.io/files/63faba13bd096f001c1de659.mov
演習課題：ホルムアルデヒドの電子励起状態の解析
ホルムアルデヒド（$ \mathrm{H_2CO}）の電子励起状態について、以下の問いに答えてみてください。WebMO（講座用）を用いて量子化学計算を実行する際には、基底関数は 6-31G(d) 、計算方法は、基底状態については密度汎関数理論（DFT）法、励起状態については時間依存密度汎関数（TD-DFT）法を用いて、汎関数は B3LYP とします。 問１
この分子の吸収スペクトルにおいて観測される電子遷移のうち、最もエネルギーの低い（波長の長い）と予想される遷移に対応する吸収波長 $ \lambda_\mathrm{abs} の値を nm 単位で答えてください
$ \lambda_\mathrm{abs} = 303 nm
問２
この分子の吸収スペクトルにおいて観測される電子遷移のうち、最もエネルギーの低い（波長の長い）と予想される遷移に対応する吸収強度（振動子強度） $ f の値を答えてください
$ f = 0
問３
この分子の吸収スペクトルにおいて観測される電子遷移のうち、最も光吸収の程度が大きい（吸収バンドが高い）と予想される遷移に対応する吸収波長 $ \lambda_\mathrm{abs} の値を nm 単位で答えてください
$ \lambda_\mathrm{abs} = 136 nm
問４
この分子の吸収スペクトルにおいて観測される電子遷移のうち、最も光吸収の程度が大きい（吸収バンドが高い）と予想される遷移に対応する吸収強度（振動子強度） $ f の値を答えてください
$ f = 0.16
質疑応答
事前にいただいたご質問
Q１．cis/trans 異性化反応の遷移状態の探索がうまくいかない。どうしたら？
異性化の「反応座標に沿った最小エネルギー経路（Minimum Energy Path; MEP）」を探索し、MEP 上で最もエネルギーが高くなる点を遷移状態探索の初期構造に用いるという方法があります
シンプルな反応過程であれば、Gaussian の Scan という方法で、MEP を探索することができます → 実践例①
Ｑ２．cis/trans 異性体の励起状態の物性研究にはどのようなものがある？
二重結合部位が光励起に伴って cis/trans 異性化する場合、消光（quenching）を伴うことがあります
この消光が分子の凝集状態によって制御される凝集誘起発光というものがあります
Ｑ３．化学反応の解析でどの範囲まで解析できるのか？
Gaussian のみを用いる解析では、始状態と終状態のエネルギー差、始状態 or 終状態と遷移状態のエネルギー差（活性化エネルギー）など、化学反応の「静的なふるまい」を求めることが多いように思います
始状態と終状態のエネルギー差からは、ボルツマン分布則に基づいて、両状態に対する「平衡定数」を算出することができます
活性化エネルギーの値からは、遷移状態理論に基づいて、対応する過程の「反応速度定数」を算出することができます
Ｑ４．遷移状態は、反応物と生成物の指定だけで解析できるのか？
遷移状態の探索で「QST2」という方法を用いると、「反応物」と「生成物」の２つを指定することで、遷移状態の探索を始めることができます → 実践例②
ただし、複雑な過程の場合はうまくいかないことが多いです
複雑な過程の場合には、「QST3」という方法を用いて、「反応物」と「生成物」に加えて、「探索をはじめる初期構造」を指定することで、収束性が良くなります
実践例①： 最小エネルギー経路（MEP）を Scan 法で決めてから、遷移状態を探索する
ブタジエンの cis-trans 異性化反応を例に
https://scrapbox.io/files/63fb4375f392ca001b16da41.mov
実践例②：QST2 法を用いて、始状態と終状態の構造だけを用いて、遷移状態を探索する
ブタジエンの cis-trans 異性化反応を例に
https://scrapbox.io/files/63fbfb38bcccfa001ce36f0c.mov
当日のご質問
さいごに
WebMO について
この講座の終了後、２週間（2023年3月14日まで）は利用していただくことができます
継続的な利用・制限なしでの利用をご希望の場合には、共同研究・技術相談などをお受けすることができます
共同研究・技術相談
計算化学（分子動力学・量子化学・機械学習）に関して、共同研究や技術相談をお受けしています
興味・関心を持っていただけましたら、山本までお問い合わせください
企業様とのこれまでの実施例
共同研究・受託研究
100 〜 200 万円 / 年
技術相談
2 〜 5 万円 / 時
事前相談は無償でお受けできますので、お気軽にお問い合わせください
オンラインにも対応いたします
大学など、学術機関同士で共同研究するときには、将来的に論文を共著で発表すること、競争的資金に応募することなどを目指して取り組むことで、お互いに win-win になればと願っております
お問い合わせ
norifumi.yamamoto (at) p.chibakoudai.jp
(at) を @ に変えてください