3E11ペプチドの振動モードの諸性質とハロゲン結合の構造的特徴の相関の解析(静岡大教育)鳥居肇

Analysis of the Relation between the Properties of the Vibrational Modesof the Peptide Group and the Structural Characteristics of Halogen Bonding(Shizuoka University) Hajime Torii
序
ハロゲン結合とは，共有結合に与るハロゲン原子が，他分子の負電荷を持つ原子と引力的に相互作用する分子間相互作用である。生体分子系では，ハロゲン原子を有するリガンドとペプチド主鎖ないし側鎖の間の相互作用系などが知られている。
このうち，ペプチド主鎖の C=O が関わる場合については，ペプチド基のアミドⅠモードの赤外強度が増大することを以前に報告した [1]。ペプチド基の代表的モデル化合物である Nメチルアセトアミド (NMA) と多様な含ハロゲン分子の会合体を対象に，アミドⅠモードの赤外強度を計算すると，孤立分子では 314.7 km mol-1［B3LYP/6-31+G(2df,p)，以下同様］であるのに対し，最大限 835.4 km mol-1（NMA...Br2 の場合）まで増大するという結果が得られる。
この増大の程度は水素結合の場合（NMA...3H2O，増大は+213.1 km mol-1）に比べて遥かに大きく，分子間電子移動の分子振動による変調が重要な役割を果たしている [1]。
この計算では，ハロゲン原子がペプチド基の O=C–N 平面内にあるものとしており，ハロゲ
ン結合形成に伴う双極子微分ベクトルの変化の向きは，X...O=C–N２面角（0º or 180 º）と良い相関を示している [1]。しかし，実際の生体分子系におけるハロゲン結合では，C–X...O は
ほぼ直線形であるものの，X...O=C–N２面角には柔軟性があり，0º や 180º から大きく外れるものが多いことが知られている [2]。例えば，cyclin-dependent protein kinase 2 (CDK2) と4,5,6,7-tetrabromobenzotriazole (TBB) の会合体 (PDBID: 1P5E) [3] では，Br...O=C–N２面角はおよそ–40 º および 130º となっている。そこで，このようなハロゲン結合の構造的特徴が，ア
ミドⅠモードの赤外強度をはじめとするペプチドの振動モードの諸性質に，どのように影響するかは，興味深いところである。
本研究では，NMA...TBB 会合体（など）を対象とした計算により，この点について理論的検討を行った。
結果と考察
ペプチド基の O 原子を中心に C=O 方向を z 軸とする極座標系（θ–φ，但し Br...O=C–N２面
角が 0º において φ = 0º）を定義し，O と直接的に相互作用する Br の θ および φ を固定してNMA...TBB 会合体の構造最適化を行い，アミドⅠモード（正確には，交換可能プロトンを重水素化した会合体のアミド I’モード）の振動数と赤外強度を計算した［B3LYP/6-31+G(2df,p)］。
結果を図１に示す。振動数は，θ が小さい領域において θ と良い相関を示す一方［図１(a)］，赤外強度は (θ, φ) = (30º, 180º) を中心とした分布をとる［図１(b)］ことがわかる。また，他の構造パラメーターとの相関を検討したところ，振動数は（予想されるとおり）C=O 結合長と良い相関を示したが［図２(a)］，赤外強度にはそのような相関は見られず（図省略），むしろO...Br 結合長との相関が良く［図２(b)］，またハロゲン結合形成に伴う双極子微分ベクトルの，変化は，θ ≤ 60º の範囲では大きさはあまり変わらず（1.055 ± 0.097 D Å-1 amu-1/2，図省略）
O...Br に近い向きをとった［図２(c)］。これらの結果は，アミドⅠモードの振動数と赤外強度から，ハロゲン結合に関わる分子どうしの相対的配置の情報が得られることを示している。
振動数が C=O 結合長および θ と良い相関を示し，φ にはあまり依存しないことは，C=O について軸対称なプロパティが振動数シフトに関与していることを示す。これは，振動数シフ
トを当該モードの双極子微分と相手分子由来の電場の相互作用によるものとして記述する理論モデルの結果と類似しているが，振動数シフトの θ 依存性が逆符号であり，詳細なメカニ
ズムについては今後の検討を要する。一方，ハロゲン結合形成に伴う双極子微分ベクトルの変化が O...Br に近い向きをとることは，これが分子間電子移動の分子振動による変調に由来するという理論解析の結果 [1] を支持する。
1735

(a)

frequency / cm

-1

1734

1730
1728

(a)

1726
1724

-1

1730

1.22351.2240

1725

0

20

40

60

80

30

60

theta / deg.
60

30

0

90

180

0

(b)

IR intensity / km mol

1720

90

1.22451.22501.2255r(C=O) / Å-1

frequency / cm

1732

600
500
400

(b)
300

600

30

150

340

phi / deg.

120

60

90

図１：NMA...TBB 会合体（交換可能プロトンを重水素化したもの）のアミド I’モードの (a) 振動数と (b) 赤外強度の Br 原子角度位置 (θ–φ)

angle between diff. dipole deriv.
and O...Br / deg.

2.90

2.95

3.00

3.05

3.10

3.15

r(O...Br) / Å
20

(c)

15
10
5
0
0

20

40
60
theta / deg.

80

依存性。後者における等高線は 20 km mol-1 おき

図２：NMA...TBB 会合体（交換可能プロトンをに描かれており，600 および 340 km mol-1 の等高重水素化したもの）のアミド I’モードの (a) 振線のみ赤外強度値を明示している。赤色の×印

動数と C=O 結合長，(b) 赤外強度と O...Br 結合

は，ハロゲン結合会合体が最適化構造として得長，(c) ハロゲン結合形成に伴う双極子微分ベク

られなかった角度位置を示す。

トルの変化の向きと θ，の関係。

[1] H. Torii, J. Chem. Phys. 133, 034504 (2010). [2] A. R. Voth, P. Khuu, K. Oishi, and P. S. Ho, Nature Chem. 1,
74 (2009). [3] E. De Moliner, N. R. Brown, and L. N. Johnson, Eur. J. Biochem. 270, 3174 (2003).

3E12
トンネル効果と反応物錯体形成が引き起こす抗酸化反応の活性化：
ユビキノール-ビタミン E 間の水素移動反応(京大院・理)○稲垣 泰一、山本 武志

Tunneling and donor-acceptor binding accelerate antioxidant reaction:hydrogen transfer reaction of ubiquinol and vitamin E.
(Kyoto University) ○Taichi Inagaki and Takeshi Yamamoto【序論】ビタミン E による生体膜中の抗酸化反応は単純な水素移動反応である。その反応によって生体は過酸化ラジカルから保護されている。水素移動によって抗酸化作用を失ったビタミ
ン E ラジカルは他から水素を受け取ることでビタミン E に再生される。水素ドナー分子の一つとしてユビキノールが挙げられるが、この水素ドナー分子によるビタミン E 再生反応は速やかに行われなければならない。なぜなら、ビタミン E ラジカルを放置すればいずれ抗酸化反応の逆反応が起こり、再び過酸化ラジカルが生成されてしまうためである。従って、ビタミン E 再生反応の反応速度は非常に重要な物理量である。
ユビキノールによるビタミン E 再生反応の反応速度は長岡らによって実験的に調べられた。彼らはユビキノール-10 (UQH2) とビタミンE 誘導体 (5,7-diisopropyl tocopheroxyl radical; Toc) (図 1)をそれぞれのモデルとして用いて溶液内実験を行った 1。その結果、非常に大きな速度論的同位体効果が得られたため、彼らは水素のトンネリン
グによって反応速度の加速がなされていると主張した。
本研究は長岡らの実験に理論的考察を加えるとともに、ビタミン E
再生反応の加速に重要な因子を明らかにすることを目的とする。長岡らの実験に則してエタノール中 UQH2-Toc 水素移動反応を対象とし、図 1. UQH2 (上)、Toc (下)反応メカニズム、トンネリング、反応速度を理論的に調べる。
【計算方法】電子状態計算 (DFT 法) と溶媒記述のための 3D-RISM 理論を組み合わせた
3D-RISM-SCF 法によって反応の自由エネルギー面を作成した。そして、その自由エネルギー面の情報をインプットとし、variational transition state theory with multi-dimensionaltunneling (VTST/MT) 法を用いて反応速度を計算した。ここでは、平衡溶媒和と断熱的電子移
動反応を仮定している。得られる速度定数は k =κkVTST の形で表現され、トンネリングは透過係数κの形で速度定数に寄与する。

【結果・考察】自由エネルギー面作成のため、まず溶質構造の主な停留点を決定した。その結果、反応中心近くに大きな置換基があるにもかかわらず、遷移状態 (TS) は非常に短いドナー
酸素-アクセプター酸素間距離 (R(Od-Oa) = 2.34 Å) を持つことがわかった。これらの構造をもとに得られた透過係数κの大きさは 4.1×103 となった。これはかなり大きな値であり、トンネ
リングによってこのモデル抗酸化反応は 4000 倍以上も加速されるということを示している。この大きなトンネリングの原因は、非常に短い R(Od-Oa)の構造を TS 付近で取り得ること、そしてそれが反応障壁の薄さ、トンネリング経路の短縮につながるためと考えられる。反応速度の評価の際には、様々な DFT 汎関数や波動関数理論を用いて反応に伴うエネルギー変化を検証した。そこでは反応物錯体が安定化されることで直接的に反応障壁が下がることが確認された。
そして、得られた速度定数 k は 3.5×104 M-1s-1 となり実験値(4.1×104 M-1s-1)と良く一致した。
以上の結果から、トンネリングと反応物錯体の安定化という 2 つの効果は、モデル抗酸化反応
はもとより、生体膜中の抗酸化反応でも反応速度上昇に大きく寄与し、重要な要素となっていることが示唆された。実験結果との他の比較としては活性化エネルギーの評価を行った。本研究では VTST/MT 法による速度定数の表式に従って解析的に活性化エネルギーを求めた。得られた活性化エネルギーは 1.0 kcal/mol 以下となり、実験値と定性的に良く合う結果となった。この値は非常に小さく、トンネリングと反応物錯体の安定化によってもたらされる。また、活性化エネルギーは

図 2. 反応プロファイルと活性化エネルギー。

平均トンネリングエネルギーとほぼ一致すること

MO は各構造の HOMO を表す。

も示された。
また、本研究では反応の電子状態と絡めた反応メカニズムの解析も行った。水素移動反応はし
ばしば、水素原子移動 (HAT) 反応とプロトン共役電子移動 (PCET) 反応に分けられるが、フェノール-フェノキシルラジカル系の水素移動反応はそのメカニズム分類で多くの注目を集めている 2。今回の解析では高精度波動関数理論を用いて基底状態、励起状態の計算を行っている。
その結果、UQH2-Toc の系では、それらが混合した反応メカニズムを持つことがわかった。HAT反応よりも電子的非断熱性の強い PCET 反応が多く混合されることでトンネリングの増大につながる可能性があることが示唆される。
[1] S. Nagaoka, et al., Chem. Phys. Lett. 287, 70-74 (1998), S. Nagaoka, et al., J. Phys. Chem. B 104,856-862 (2000), A. Ouchi, et al., J. Phys. Chem. B 114, 6601-6607 (2010).
[2] J. M. Mayer, et al., J. Am. Chem. Soc. 124, 11142-11147 (2002), J. H. Skone, et al., J. Am. Chem.

Soc. 128, 16655-16663 (2006), A. Cembran, et al., J. Chem. Theory Comput. 8, 4347-4358 (2012).
[3] T. Inagaki and T. Yamamoto (submitted).

3E13

3D-RISM-SCF 法によるメロシアニンの電子状態に対する溶媒効果の解析(九大院理) ○田中佑一, 吉田紀生, 中野晴之

3D-RISM-SCF analysis of the solvent effecton the electronic structures of merocyanines(Kyushu Univ.) ○Yuichi Tanaka, Norio Yoshida, and Haruyuki Nakano【序論】
ソルバトクロミズムとは、溶媒の種類によって溶液の色が変化する現象のことである。この現
象は溶媒の種類によって溶質の電子状態が変化し、それに伴って励起エネルギーがシフトすることによって生じる。このときのシフトの大きさは大まかには溶媒の極性によって説明されるが、例えばメタノールとアセトニトリルのように誘電率が比較的近い溶媒でも異なる色を示すなど、溶媒の極性だけでは十分に説明できない場合もある。そのため、溶質-溶媒間相互作用を分子論的に調査することは極めて重要である。
そこで本研究では、直鎖状のモデルメロシアニンである streptopolymethinemerocyanine(SPMC, 図 1) の溶液中での励起エネルギーを three-dimensional reference interaction sitemodel self-consistent field (3D-RISM-SCF) 法および時間依存
密度汎関数理論 (TD-DFT) によって求め、励起エネルギーの溶媒依存性と鎖長依存性について系統的に調査した。また、その要因について溶媒和構造から議論した。
なお、
これまでも RISM-SCF
法を用いた調査を行ってきたが[1]、本研究では 3D-RISM-SCF
法を用いて、より詳細な溶媒和構造の解析を行った。

図 1 SPMC【計算方法】SPMC (n=1-8) を溶質として水、メタノール、アセトニトリル溶媒中での計算を行った。電子状態の手法として密度汎関数理論 (DFT) を用い、基底関数として aug-cc-pVDZ を、汎関数として LC-BOP を使用した。SPMC を孤立状態で構造最適化した分子構造を用いて、3D-RISM-SCF
計算および TD-DFT 計算を行い、溶液中での垂直励起エネルギー、溶媒和構造等を求めた。

【結果と考察】図 2 に SPMC の垂直励起エネルギーを示した。まず溶媒依存性に着目する。鎖長が短いとき (例えば n=1 のとき) は–*励起エネルギーはレッドシフトし、n–*励起エネルギーはブルーシフトする。また、シフトの大きさは水中>メタノール中>アセトニトリル中の順となる。一方、鎖長が長いとき (例えば n=8 のとき) は、
水中とメタノール中では鎖長が短いときと同様に–*励起エネ
ルギーはレッドシフトし、n–*
励起エネルギーはブルーシフト
するが、アセトニトリル中では–*励起エネルギーはブルーシフ
トし、n–*励起エネルギーは気相中からほとんどシフトしない。
次に鎖長依存性に着目する。

図 2 SPMC の垂直励起エネルギー気相中でも溶液中でも鎖長を伸ばすにつれて–*・n–*励起エネルギーともに減少する。しかし、その減少の度合いは溶媒の種類に依存して異なっている。その結果、–*励起エネルギーと n–*
励起エネルギーの大小が気相中では n=3 と 4 の間で入れ替わっていたのが、水中では n=1 と 2の間、メタノール中では n=2 と 3 の間で入れ替わるようになる (アセトニトリル中では気相中と同様に n=3 と 4 の間で入れ替わる) 。
よって、これらの結果は SPMC の励起エネルギーが溶媒と鎖長の両方に依存することを示している。
このような依存性を示す要因を溶媒和構造により考察する。図 3 に水中での SPMC (n=1) の周りの溶媒和
構造を、図 4 にアセトニトリル中での SPMC (n=1) の周りの溶媒和構造を示した。これらの図より、水中では SPMC は両端 (カルボニル基の酸素原子とアミノ基の水素原子) で溶媒と水素結合を形成できるが、アセトニトリル中では SPMC はアミノ基側のみで溶媒と水素

図 3 SPMC (n=1) の水中での溶媒和構造結合を形成していることがわかる。そのため、溶質-溶媒間相互作用は水中の方がアセトニトリル中よりも大きくなり、その結果、励起エネルギーの溶媒和シフトの大きさが水中>アセトニトリル中となることが説明される。
鎖長が長いときの溶媒和構造による考察については発表当日に行う。
図 4 SPMC (n=1) の【参考文献】

アセトニトリル中での溶媒和構造[1] Y. Tanaka, N. Yoshida, H. Nakano, Chem. Phys. Lett. in press.

3E14
線形応答近 似と新規静電場フ ィッティング法を組 み合わせた LR-RISM-SEDD 法の開発 	
 (名大 WPI-ITbM)	
 横川大輔	 Development of LR-RISM-SEDD by combining linear
re sponse approximation and new fitting approach ofelectrostatic pote ntial(WPI-ITbM) Daisuke Yokogawa【 序 論 】近年の計算機の発達により、周りの環境を考慮した、より実在に近い系での量子化学計算が行われるようになってきた。特に溶液内での量子化学計算は、化学反応、光化学、生物学など様々な領域で必要とされている。これまでに連続誘電体モデルや分子動力学法と量子化学計算を組み合わせた QM/MM-MD 法など様々な手法が開発されているが、精度と計算効率を両立させる困難さから、未だに多くの研究者によって新たな理論開発が行われている。
	  我々は積分方程式理論の一つである Reference Interaction Site Model (RISM)法と量子化学計算を組み合わせた RISM-SCF-SEDD 法の開発を行ってきた。本研究では計算精度と計
算 コ ス ト の 大 幅 な 向 上 を 目 指 し 、 RISM-SCF-SEDD 法 に 線 形 近 似 を 適 用 し た
LR-RISM-SEDD 法の開発を行った。
溶液内量子化学計算で常に問題になることは、着目する溶質周りに存在する無数の溶媒分子を如何に取り扱うかである。RISM-SCF-SEDD 法[1]では、反応中心を量子化学計算で解析を行い、溶媒分子は統計力学的手法の一つである Reference Interaction Site Model (RISM法)を用いて取り扱っていた。この手法は溶液内の様々な化学反応、光化学に適用されてはいるものの、計算コストが高い、分極が大きな分子で計算が不安定になる等の問題点があった。
そこで、本研究では線形応答近似(LR)と新規静電場フィッティング法を導入することで、この両者の問題点の克服を試みた(LR-RISM-SEDD 法[2])。
【 方 法 】 本手法では電子密度ρを、補助基底χを用いて次式で近似する。

d は展開係数であり、次式を最小化することで決定した。

ここで、Ne は電子数、θは Heaviside 関数である。μは近似した電子密度が負の値をとらないように、大きな値(本研究では 500)を用いている。すべての計算は本手法を新たに組み込ん

だ GAMESS プログラムをTable	
   1.	
   Computational	
   times	
   required	
   in	
   LR-­‐RISM-­‐SEDD	
   and	  用いて行った。

RISM-­‐SCF-­‐SEDD	
  (unit:	
  second).	
  【 結 果 と 考 察 】本手法の計

	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   LR-­‐RISM-­‐SEDD	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   	
   RISM-­‐SCF-­‐SEDD

算コストについて検討する

#	
  of	
  cores	
   	
   	
   	
   1	
   	
  

6	
   	
  

12	
   	
   	
  

1	
   	
  

6	
   	
  

12

た め に H2O 、 H2CO3 、

H O	
   	
  

32	
   	
  

22	
   	
  

22	
   	
  

193	
   	
  

193	
  

	
   194p-nitroaniline (PNA)につい

H CO 	
   	
  

241	
   	
   189	
   	
  

186	
   	
  

1880	
   	
  

1967	
   	
  

2003

て RISM-SCF-SEDD 法と本

PNA	
   	
  

2538	
   	
   1895	
   	
   1904	
   	
  

2
2

3

28194	   	
   28480	
   	
   29049

手法で、計算コアを変えつつ
計 算 時 間 の 測 定 を 行 っ た (Table 1) 。 1 コ ア を 使 用 し た 場 合 、 本 手 法 の 計 算 時 間 は
RISM-SCF-SEDD 法よりも約 10 分の 1 であることがわかった。この計算時間の差は用いる
コア数の増加とともに大きくなり、PNA を 12 コアで計算した場合、本手法の計算時間はRISM-SCF-SEDD 法の約 15 分の 1 となることがわかった。
計算精度について、Fig. 1 に示す 2-phenacylpyridine の ketoimine (KH)、enaminone (MH)、enol (EH)の相対安定性を計算することで検討を行った。構造最適化、振動数計算は B3LYP/6-31G*、
エネルギー計算は 6-311+G**を用いて行った。
CCSD(T) 法 に つ い て は 6-31+G* で 計 算 を 行 い
MP2/6-311+G**で得られたエネルギーで補正を行った。

Fig. 1. Tautomers of 2-phenacylpyridine.

	  Table 2 に、KH、MH、EH の相対安定性を気相中、水溶液中で計算した結果を示した。
水溶液中の計算では本手法による結果だけでなく、PCM 法による結果も示した。
	
  CCSD(T) 法 と本手法を組み合わせることで、実験
値を良く再現できることがわかった。
発表では、構造変化や溶媒構造変化を示しながら、詳

Table	
   2.	
   Relative	   free	
   energies	
   computed	
   in	
   gas	   and	   aqueous	
   phases	
   (unit:	
  kcal/mol).	
   	
  
	
   	
   	
   	
   	
   	
   	
   gas	  phase	
   	
   	
   	
   	
   	
   	
   	
   aqueous	
  (PCM)	
   	
   	
   aqueous	
  (this	
  work)	
   	
   	
   	
   exp.[3]	
  
	
   	
  

B3LYP	   	
   CCSD(T)	
   	
  B3LYP	   	
   MP2	   	
  B3LYP	   	
   CCSD(T)	
  

KH	
   	
  

0.0	
   	
  

0.0	
   	
  

0.0	
   	
  

0.0	
   	
  

0.0	
   	
  

0.0	
   	
  

0.0	
  

MH	
   	
  

-­‐0.2	
   	
  

8.7	
   	
  

-­‐0.6	
   	
  

8.6	
   	
  

-­‐5.5	
   	
  

1.3	
   	
  

1.4	
  

EH	
   	
  

-­‐2.7	
   	
  

2.6	
   	
  

0.0	
   	
  

4.2	
   	
  

1.1	
   	
  

3.8	
   	
  

2.7	
  

細について議論する予定である。
【参考文献】
1. D. Yokogawa, H. Sato, and S. Sakaki, J. Chem. Phys. 126, 244504 (2007).
2. D. Yokogawa, J. Chem. Phys. 138, 164109 (2013).
3. A.R.E. Carey, et al., J. Chem. Soc., Perkin Trans. 2, 2285 (1993).

3E15
フラグメント法を用いた化学反応の理論考察（産総研・ナノシステム 1，東工大・生命理工 2，神戸大院システム情報 3）○Fedorov
Dmitri1，中田 浩弥 2，北浦 和夫 3
Theoretical study of chemical reactions with fragment-based methods(NRI, AIST1; Tokyo Inst. Techn.2; Kobe U.3) ○Dmitri G. Fedorov1, Hiroya Nakata2,
Kazuo Kitaura3【序】
量子化学計算法は汎用性乃至系統的改良性を持てど、計算量は系の規模に対して、三乗以上程大きくなる為、巨大系への応用は難しい。分子系をフラグメントに分割する方法、即ちフラグメン
ト法[1]で計算は部分的に行われる為、規模の依存を望ましくば一乗まで削減出来る。フラグメン
ト法として、フラグメント分子軌道(FMO)法[2,3,4]を中心に、有効フラグメント分子軌道(EFMO)
法[5]と有効フラグメントポテンシャル(EFP)の QM/EFP[6]を化学反応に適用した。
【方法】
1. FMOFMO 法では全系のエネルギーを下記の式でフラグメントのエネルギーから得られる。
N

N

I =1

I >J

E = ∑ E I + ∑ (E IJ − E I − E J )EI と EIJ はフラグメント I とフラグメント対 IJ のエネルギーである。
解析二次微分[3]を用いて、虚数振動数を追跡し、遷移状態を探索する。現在、FMO 法の解析二次
微分は静電場近似を使わない FMO2-RHF のためにしか開発されず、計算量が大きいが、必要なメモ
リは系の規模の二乗となるので、巨大系でも適用は可能である。
化学反応を解析する為に、制限開殻対相互作用解析法(RO-PIEDA)を開発した。その方法では、フラグメント間の相互作用を四つの成分に分割する。即ち、静電(ES)、交換反発(EX)、電荷移動と高次項(CT+mix)と分散力(DI)の成分を以って、化学反応による相互作用の変化を見て、反応を解明する。
2. EFMOEFMO は FMO によく似た計算法で、分極の扱いとその電子状態への影響は主に異なる所である。
EFMOの解析二次微分は開発されていない為、遷移状態を探索する為に、別な作戦を執った。固定領域(FD)の EFMO/FD 法を開発し、それで化学反応中心の原子の位置を最適化した。最適化中各分子構
造の電子状態を化学反応中心含有可分極領域に限り解く為、計算量を大分削減する。遷移状態を探索する為に、予め反応経路を決めて、その経路に沿って反応座標を固定しながら、他の自由度を最適化する（部分的拘束構造最適化）。
3. QM/EFPFMO と EFMO の実用的動力学模擬の計算量が大きい為、より高速な計算法は当面必要である。水分子を有効ポテンシャル(EFP1)で扱う QM/EFP を用いて、複製交換動力学(REMD)の模擬を実装した。
高速に並列計算を行うに、全 CPU コアを班に分ける二段階並列化(GDDI)法と組み合わせ、各複製計算を専用班で行って、高い並列化効率の模擬を実現した。

【結果】全ての計算では GAMESS と数十コアの PC cluster を用いた。
1. FMO の結果
二次微分を用いて、CH3X+OH––→CH3OH+X– (X=F, Cl, Br)の SN2 反応の遷移状態を FMO-RHF と分割を使わない RHF で行い、遷移状態構造や反応障壁を比較した。

(a)FMO-RHF で得た虚数振動(b)RHF の虚数振動

別な課題として、(ALA)3 の水素抽出反応に RO-PIEDA を適用し、ラジカル安定化に伴い、抽出される水素の位置による影響を解析した。
2. EFMO の結果
FD の近似法(FMO/FDD と EFMO/FDD)を用いて、酵素反応の経路を算出した。EFMO 法と FMO 法の結果はほぼ一致した。幾つかの方法で得られた反応障壁を比較して、RHF は大きく違う一方、MP2 とDFT の結果はよく一致すると分かった。

Potential of mean force(kcal/mol)

3. QM/EFP の結果
過酸化水素分子を 290 水分子に浸し、

9
HF/3-21G//TIP5P

8

水和中異性体化反応の動力学の模擬とHF/6-31G*//TIP5P
MP2/6-31G*//TIP5P

7

して、1ns の REMD と HOOH 二面角度をMP2/aug-cc-pVDZ//TIP5P

6

MP2/aug-cc-pVDZ

反応座標とした傘型抽出法(umbrella

1ns REMD-32 Replicas

5

2ns REMD-32 Replicas

sampling)の 100 ps の MD/US を比較し

1ns REMD-64 Replicas

43た。真空中の結果（黒点線）は水和中

2と大分違う一方、REMD と MD/US の反応

1

経路はよく一致している。trans 型の

0
0

30

60

90

120

150

Torsional Angle, τHOOH (Degree)180

構造は双極子モメントを持ち、水和中安定化される。

【結論】
計算量や汎用性或いは精度が異なるフラグメント法は化学反応追跡と動力学模擬に有用であると明らかにした。現象と計算機性能に合わせて方法を選択すれば、現実的応用が可能である。
参照[1] M. S. Gordon, D. G. Fedorov, S. R. Pruitt, L. V. Slipchenko, Chem. Rev. 112 (2012) 632.
[2] http://staff.aist.go.jp/d.g.fedorov/fmo/main.html[3] H.Nakata,T.Nagata, D.G.Fedorov, S.Yokojima,K.Kitaura,S.Nakamura, J.Chem.Phys.138(2013) 164103.
[4] M.C. Green, D.G. Fedorov, K.Kitaura, J.S. Francisco, L.V. Slipchenko, J.Chem.Phys.138(2013) 074111.
[5] C. Steinmann, D. G. Fedorov, J. H. Jensen, PLoS ONE 8 (2013) e60602.
[6] D. G. Fedorov, Y. Sugita, C. H. Choi, J. Phys. Chem. B 117 (2013) 7996.

3E17

NO3 分子の振動回転スペクトルの解釈に関する
計算化学からのアプローチ（お茶大理 1，産総研 2）○平野 恒夫 1,2, 長嶋 雲兵 2
An approach from computational chemistry to the interpretation of the rotationand vibration spectra of NO3 molecule(Ochanomizu U. 1, AIST 2) ○Tsuneo Hirano1,2, Umpei Nagashima2NO3 の振動スペクトルの帰属に関して廣田グループと Stanton グループで 20 年以上にも及ぶ論争が続いている。最大の問題点は、1492 cm-1 の band の帰属をν3 と考えるか 1)、ν1 + ν4 2) または
ν3 + ν4 3)と考えるか、の違いである。従来、Stanton らの CCSD 系の計算では D3h 対称の構造がグ
ローバルミニマムになることと、回転スペクトルの解析結果は D3h 対称の構造を示唆していることで、D3h 構造を仮定して振動・回転スペクトルを統一的に解釈しようとなされてきた。また、平
衡構造近辺のポテンシャル曲面を計算して D3h 対称の構造の周辺にある幾つかの local minima を
dynamic に動き廻っているとするモデルをたてて解釈しようとしたグループ 4) もある。 我々は、励起状態からの摂動があって、その上 dynamic な構造変化をしているのであれば、multi-reference
framework に準じた計算が必要と考え、MR-SDCI+Q の計算を行ったところ、D3h 構造は 2 次の遷移状態で、グローバルミニマムはそれより約 800 cm-1 低い C2v 構造である
ことが分かった。
まず、RCCSD(T)や DFT など単一配置（Single reference framework）の計算を試みたところ確かに D3h 対称の構造が global minimum となっていた。
つぎに、多配置の計算として、CASSCF とその波動関数に基づく多配置の
CI、つまり MR-SDCI＋Q（Q は Davidson correction）を行ってみた。用いた基底関数は aug_cc_pVQZ である。計算には MOLPRO 2010.1 を使った。

Fig. 1 C2v Structure

NO3 の構造式は形式的に Fig. 1 のように書けるので、とりあえず C2v 対称の構造になるように
z-matrix で指定して、nosym 条件の下で最低根を計算した。結果を Table 1 に示す。
D3h 構造を保つようにして構造最適化した構造（Table 1

Table 1 Stationary structures from ab initio calculationsの D3h ）は、CASSCF でもSymmetry

MR-SDCI+Q ででも、2 個の

CASSCF (1s closed)

虚の振動数を有する 2 次のD3h

1.24251.2425120.0

600

遷移状態（TS）であった。

C2v

1.34571.2095114.6

0C2v の構造を指定して求めた最適化構造には(a)と(b) 2 個

D3h

典的な構造式（Fig. 1）に相C2v

C2v 構造の(a)と(b)をつなぐr(NO1)

theta

Relative E /cm-1

MR-SDCI+Q (1s core, valence closed)が存在した。その内、(a)は古当する構造になっている。

r(NO)1.2185

1.2185

120.0

997

（a)

1.3242

1.1915

114.4

0

TS

1.21901.2185

120.0

842

(b)

1.18281.2466126.1

178

minimum energy path (MEP)
を Fig. 1 の theta を変数として求めたところ、Fig. 2 のようになった。MEP に沿って r(NO)と r(NO1)が入れ替わるように、大きく変化し、両者が交差するところが Table 1 の C2v TS に相当する。

この TS は、(a)からは 842
cm-1, (b)からは 664 cm-1 の高さなので、もし、このポテ
ンシャルエネルギー曲面の計算が正しければ、振動ス
ペクトルにはトンネル効果が見える筈である。
この計算では、任意の 1
本の N-O 結合を C2v の主軸になるように選んで nosymの計算をやっている。分子全体では三本の N-O があるので、dynamic に構造を変化させつつ振動・回転がFig. 2 Potential energy E at the MR-SDCI+Q level and N-O bondlength r(N-O)s along the MEP. See Fig. 1 for the definitions oftheta, r(NO) and r(NO1).

励起されるプロセスを考えると見かけ上極めて複雑なポテンシャルエネルギー曲面上で振動・回転が起こっていることになる。
予想外に複雑なポテンシャル曲面に驚いているが、こうなると、仮に面外振動は無視するとしても、「C2v の制約の下での分子の変形」という条件を緩めて Cs での構造を探る努力をする必要があると考えて計算を続けている。振動・回転エネルギーの議論に入る前にまず、しっかりした
ポテンシャルエネルギー曲面を構築することが緊要の課題であると考える。
測定のタイムスケールを考えると、振動スペクトルは C2v 構造に対応したものに、回転スペク
トルは時間平均されて D3h 構造に対応するものになっているものと考えられる。現時点の計算結果はまだ予備段階程度のものであるが、積年の論争の解決の糸口になるのではないかと思って報告した。

例えば、E. Hirota, T. Ishiwata, K. Kawaguchi, M. Fujitake, et. al., J. Chem. Phys., 107, 2829 (1997) .
J.F. Stanton, J. Chem. Phys. 126, 134309. (2007).
3
K. Kawaguchi, N. Shimizu, R. Fujimori, J. Tang, T. Ishiwata, I. Tanaka, J. Mol. Spectrosc., 268, 85 (2011).
4
M. Mayer, L.S. Cederbaum, H, Köppel, J. Chem. Phys., 100, 899 (1994)
１

2

3E18
ZnOH の計算分子分光学: FeOH, CsOH との比較

（産総研 1，Wuppertal U. 2，お茶大理 3）○長嶋 雲兵 1，Per Jensen2, 平野 恒夫 1, 3
Computational Molecular Spectroscopy of ZnOH: Comparison with FeOH and CsOH
(AIST1, Wuppertal U.2, Ochanomizu U.3) ○Umpei Nagashima1, Per Jensen2, Tsuneo Hirano1,3
我々は、既に電子基底状態の 1Σ+ CsOH と 6A’ FeOH について、それぞれ RCCSD(T)_DK3 およびMR-SDCI+Q_DK3 レベルの計算を行い、その平衡構造が CsOH では linear1)、
FeOH では quasi-linear2)であることを明らかにした。
ここで、
quasi-linear とは変角振動のポテンシャルエネルギー曲面（PES）において、直線型の構造の位置に変角振動のゼロ点振動程度の障壁を持つ分子のことをいう。波動関数の解析から CsOH の Cs-O 結合はイオン結合で、1) FeOH の Fe-O 結合もイオン結合ではあるが若干の共有結合性が認められること 2) が分かっており、さらに共有結合性が強まっていることが期待される ZnOH について、MR-SDCI+Q_DK3/[QZP ANO-RCC (Zn, O, H)]レベルの計算を行って PES を定めた後、その PES に基づいて MORBID プログラムにより振動・回転のシュレーデ
ィンガー方程式を解いて、振動・回転に関する知見を求め、未だに実験報告のない振動スペクト
ルのシミュレーションを行った。MR-SDCI+Q_DK3 の計算は MOLPRO 2010.1 を使用し、active
space は Zn の 4s、4p と O の 2s、2p、H の 1s とし、Zn の 3d からの動的電子相関も取り込んだ。
ZnOH の電子基底状態は 2A’（直線構造では 2Σ+）で、その平衡構造は bent 型で、re(Zn-O) = 1.8028
Å、re(Zn-O) = 0.9606 Å、∠e(Zn-O-H) = 114.9°であった。Fig. 1 に変角振動における minimum
energy path (MEP)を示す。 比較のために、1Σ+CsOH と 6A’ FeOH の場合の PES1,2)も描き込んで
ある。

1 +

Σ

CsOH が linear 、 6A’ FeOH がquasi-linear、2A’ ZnOH が bent 型であることが分かる。また、その変角振動の 2 次の力の定数から、この三つの分子は、それぞれ大振幅変角振動をしていることが分かった。
電子状態に関して、2A’ ZnOH の平衡構造において、89%の重みを持つ配置状態関数（CSF）は、[Ar core] (3d)10 (O2s)2 (OHσ)2 (ZnOσ)2 (O2pπ)2(ZnOσ*)1と書けることが分かった。自然軌道(NO)で現した OHσ、ZnOσ、ZnOσ* 軌道を等値面で表示すると、次頁の Fig. 2 のようになる。 OHσ 軌道
は、水の OHσ軌道とほぼ同じなので、ZnOH のOH 結合は通常の共有結合である。ZnOσ 軌道は若干の Zn の 4s 軌道と主成分のΟの 2pσから成Fig. 1. Bending Potentials along the MEP of 1Σ+ CsOH,1)
6
A’ FeOH,2) and 2A’ ZnOH; 𝜌̅ = (180° − bond-angle).

る結合性のσ軌道であり、
ZnOσ* は Zn の背後に大きく張り出した反結合性の 4s-4pσ軌道であって、両者の占有数はそれぞれ 1.96 と 1.00 なので、Zn-O 結合は、本質的にはイオン結合であるが、かなりの共有結合性をもったイオン結合であると結論できる。 共有結合は方向性を持つので、予想通り、金属-酸素イオン結合における共有結合性の多寡が linear (CsOH）か、quasi-linear （FeOH）
か、bent （ZnOH）か を決めていることが分かったことになる。また、ZnOσ* 軌道の形から不対電子は ほぼ完全に Zn 上に局在していることも分かる。
Zn および Zn を含む分子は両性（amphoteric）で,酸としても塩基としても振る舞うことが知られている。ZnOH の断熱的解離極限を調べると、
ZnOH → ZnO + H
ZnOH → Zn + OHとなり、それぞれの解離エネルギーは、42600 cm-1
(5.28 eV)、44600 cm-1 (5.53 eV)であり、確かに両性分子であることが分かった。
電子基底状態 6A’ ZnOH の平衡構造における垂直励起エネルギーは、下から 29200 cm-1、35540 cm-1
で、この二つは平衡構造が直線型である 2Πに由来する A’と A’’に相当するが、これらの電子励起状態

Fig. 2. MR-SDCI_DK3 Natural Orbitals.
Occupation numbers are in parentheses,
and NO orbital-energies (see ref. 3 for thedefinition) are under the spin-configuration.

からの振動・回転-基底状態への摂動は垂直励起エネルギーが極めて大きいので無視してよいことが分かった。
は r0(Zn-O) = 1.8078 Å、
r0(Zn-O)
MORBID の解析で、ZnOH の振動・回転の基底状態の構造（r0 構造）= 0.9778 Å、∠0(Zn-O-H) = 117°であり、この r0 構造から計算した回転定数は実験値 4)と 0.3%の誤差で一致した。ZnOH の振動スペクトルの測定は未だなされていない。OH 伸縮振動、Zn-O-H 変
角振動、Zn-O 伸縮振動の調和振動数は、それぞれ 3855、695、627 cm-1、また非調和項を含んだ
振動数はそれぞれ 3677、665、616 cm-1 と計算された。
一般に、振動と回転を分離したスペクトルとして測定することが出来ないので、変角振動の波
動関数は 2 次元以上になる。 これを変角に関する 1 次元の座標で表現すると各点は回転に関して平均化された値になるので、調和振動近似の場合ででも波動関数は原点に極大値を持たない。すなわち、平衡構造が直線型であっても振動・回転の平均構造は直線からズレた曲がった構造として測定されることになる。 例えば、実験に対応する r0 構造を計算すると、CsOH、FeOH、ZnOHの場合の直線構造からのズレは、それぞれ 17、39、63 度となる。従って、実験では実測スペクト
ルから、linear か quasi-linear か（または bent か）を区別することは容易ではない。しかし、スペ
クトル解析の初期の段階から試行計算の出来る Yamada-Winnewisser quasi-linearity parameter5) γ0 を使えば、linear では−1、quasi-linear では 0、bent では 1 に近い数値になる筈なので、その区別ができる。計算では、最初にポテンシャルエネルギー曲面が求まるので（Fig. 1）、その区別は自明で確かに linear な CsOH で−1.07、
あるが、MORBID の解析で得られた振動数からγ0 を計算してみると、
quasi-linear な FeOH で 0.10、bent な ZnOH で 0.84 となっており予想通りであった。
T. Hirano, V. Derpmann, U. Nagashima, P. Jensen, J. Mol. Spectrosc., 263, 150 (2010).
T. Hirano, U. Nagashima, G. Winnewisser, P. Jensen, J. Chem. Phys., 132, 094303 (2010).
3
MOLPRO Manual および T. Hirano, M. Amano, Y. Mitsui, et al., J. Mol. Spectrosc., 243, 267 (2007).
4
L. N. Zack, M. Sun, M. Bucchino, D. J. Clouthier, L. M. Ziurys, J. Phys. Chem. A, 116, 1542 (2012).
5
K. Yamada, M. Winnewisser, Z. Naturforsch., 31a, 139 (1976).
１

2

3E19
DC-DFTB 理論の理論展開とハイブリッド並列化
1

2

3

5

（分子研 ，早大理工研 ，早大先進理工 ，JST-CREST4 ，名大院理 ）
1,2
2,3,4
西澤宏晃 ，中井浩巳，Stephan Irle5

Theoretical extension and hybrid parallelization of DC-DFTB theory（IMS1 , Waseda Univ. RISE2 , Waseda Univ.3 , JST-CREST4 , Nagoya Univ.5 ）Hiroaki Nishizawa1,2 , Hiromi Nakai2,3,4 , Stephan Irle5
【緒言】
70.0
O(N1.09)
密度汎関数強束縛 (DFTB) 法 [1, 2] は，密度汎Total CPU time
60.0
Construct H0 and S
関数理論 (DFT) に基づく半経験的手法である。2Construct γ
50.0
Diagonalize
原子間の距離とそれらの種類によって変動するパ.]ces
O(N1.00)
Determine the Fermi level[ 40.0
ラメータ [1] を用いており，AM1，PM3 などほか
Construct densityemti 30.0の半経験的手法に比べて高精度な結果を得ることUPCができる。また，積分を行うことなくすべての要素
20.0
O(N1.52)を計算することが可能なため，その計算コストは比10.0O(N1.15)
較的低い。これまでの研究では，DFTB 法を大規O(N1.54)
0.0
0
5000
10000
15000
20000
模分子へ適用する際にボトルネックとなる対角化
Number of atomsの計算コストを削減するために，中井，小林らによFig. 1 Computatinal time for energy
り開発されてきた分割統治 (DC) 法 [3] を DFTB
calculation of C2n H2n+2 obtained by DC法に適用し，DC-DFTB 法の開発を行った。その
DFTB.
結果，バッファ領域を充分にとることにより精度を損なうことなく，対角化の計算コストが系の大きさ N に対して O(N ) となった (Fig. 1)。また，全体の計算時間も 1 次元 20,000 原子で 70 秒と大幅な計算時間の削減に成功した。
しかし，実際の大規模分子では 3 次元的な拡がりを持つために，分割された領域の数，領域に含まれる原子数が，ともに系の大きさの 3 乗で増加する。増加した対角化のコストを領域の数だけ行う必要があるため，シングルコアでの分子動力学計算は現実的でない。そこで，本研究では次世代スー
パーコンピュータ「京」での実行を踏まえて，DC-DFTB 法に対する効率的な並列化手法を開発したので報告する。また，大規模分子に対する理論的研究のために DC-DFTB 理論に基づく Hessian の実装，分子動力学計算への適用なども行ったので，併せて報告する。
【理論】
DFTB 法のエネルギーは次のように表される。

E0 =

∑

Dµν (hµν + Fµν ) + Erep =

µν

∑

1 ∑
γAB ∆qA ∆qB + Erep
2
atom

0
Dµν Hµν
+

µν

(1)

A,B

0
ここで，Hµν，γAB ，Erep は上述したパラメータにより求められる。Dµν は電子密度，∆qA は原子
A の誘起電荷である。DC-DFTB 法では全系を重なりのない部分系 α に分割し，周囲の効果を取り込むためにバッファ領域を加えて局在化領域とする。局在化領域に対して DFTB 方程式
α α α
Fα Cα
i = ϵi S Ci

(2)

α
α
α
を解くことで，部分系 α の軌道係数 Cαi を得る。F ，ϵi ，S はそれぞれ部分系 α の Fock 行列，軌道エネルギー，重なり行列である。電子数保存の制約条件から共通の Fermi 準位を決定し，密度行列を次のように構築する。
DC
Dµν
=

∑
α

α
pα
µν Dµν =

∑
α

pα
µν

∑
i

α α
fβ (ϵF − ϵi )Cµi
Cνi

(3)

ここで pαµν は分割行列，ϵF は Fermi 準位，fβ (x) は Fermi 関数である。∆qA が変化しなくなるまで式 (2) の方程式を繰り返し解き，エネルギーを求める。その後，得られた密度行列を用いてエネル
ギー勾配を計算する。エネルギー勾配計算には Yang，Lee の手法 [4] を用いた。さらに，DC-DFTB
法の Hessian への展開も行った。Hessian 表式中に現れる密度行列の微分 ∂D/∂r を，密度行列 D の導出と同様に，

]
[ α∗
DC
α
α
∑
∑
∑
∂Cµi α
∂Dµν
∂Dµν
α
α
α
α∗ ∂Cνi=pµν=pµν
fβ (ϵF − ϵi )Cνi + Cµi
∂r
∂r
∂r
∂r
α
α
i
と定義する。∂Cα /∂r は

α
∑
∂Cµi
α
α
Umi
=
Cµm
∂r
m

(4)

(5)

を満たし，Uα は部分系 α に対する coupled-perturbed Hartee-Fock (CPHF) 方程式を解くことで得られる。
【結果と考察】
本稿では DC-DFTB 法に基づく Hessian について示す。適用した系はポリエチレン (C2n H2n+2 )分子であり，計算はシングルコアで行った。DC-DFTB 法においてはユニットを C2 H2 とし，部分系
を 1 ユニット，バッファ領域を左右 8 ユニットとした。また，CPHF の収束判定は max(∂Dµν /∂b) <1.0 × 10−5 とし，収束には修正 Broyden 法を用いた。
Table 1 に n を変化させた時の平均の収束回数をTable 1 Number of iterations for CPHF示す。収束回数の最大値を 500 回としたが，通常
DFTB
DC-DFTBの DFTB 法では CPHF の収束性が悪く，n = 25の場合を除いて 500 回の時点でも 10−3 以上の変化
C50 H52
226.4330.32があった。それに対して，DC-DFTB 法ではバッ
C100 H102 Not conv.
43.70
ファ領域のみの変化を考慮しているために，収束
C150 H152 Not conv.
56.02
が非常に良い。また n = 25 のとき，通常法とのC200 H202 Not conv.
70.13
平均絶対誤差も 10−4 と精度良く求めることに成功した。
50.0
O(N1.69)
次に DC-DFTB 法の計算時間について Fig. 2
Average time per CPHF (TCPHF)
Average time per iteration (Titer)に示す。全体の計算時間を Ttot とし，1 つの座]. 40.0ce
標 に 関 し て CPHF を 解 く 平 均 時 間（TCPHF =s[
Ttot /2N ）と CPHF1 サイクルの平均時間（Titer = e 30.0mi
TCPHF /Niter ）を 示 し た 。図 に 示 し た よ う に ， tU 20.0
P
CPHF1 サイクルの計算時間は O(N 1.09 ) であり， C
10.0
計算時間も非常に小さい。現在の実装では，収束O(N1.09)
回数が系のサイズに対して線形に増加しているた
0.0
0
100
200
300
400
500
め，CPHF を解く時間は O(N 1.69 ) となっている。
Number of atoms
DIIS などの実装により収束性を改善することで，Fig. 2 Average time of CPHF and its 1
より高速な計算が可能であると考えられる。
cycle procedure obtained by DC-DFTB.
当日は Hessian 行列から得られる振動スペクト
ルや，「京」を用いた並列化性能とその並列化手法
に関しても報告する。
[1] D. Porezag, Th. Frauenheim, Th. Köhler, Phys. Rev. B, 51 (1995) 12947.
[2] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai, G.
Seifert, Phys. Rev. B, 58 (1998) 7260.
[3] M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics:
Methods and Applications (2011, Springer), pp. 97-127.
[4] W. Yang, T-S. Lee, J. Chem. Phys., 103 (1995) 5674.

3E20
化学的理解のための電子状態理論（早大先進理工、早大理工研、JST-CREST、京大 ESICB）○中井 浩巳

Electronic State Theory for Chemical Understanding(ASE, Waseda Univ.; RISE, Waseda Univ.; JST-CREST; ESICB, Kyoto Univ.) ○Hiromi Nakai
量子力学に基づいて化学を理論的な立場から研究しようとする学問、すなわち理論化学が誕生したのは 1927 年である。80 年以上の歴史において、理論化学は化学分野全体に様々な役割を果たしてきた。今後もその役割が拡大することは疑う余地もない。本講演者は、量子化学を中心とした理論化学の歴史を、1927 年から 1940 年代の第 1 期（ただし、暗黒の時代と呼ばれる 1940 年
代の発展は限定的）、1950 年代から 1980 年代の第 2 期、そして、1990 年代から 2000 年代の第 3
期の 3 つに分類して捉えている[1,2]。第 1 期は理論化学の創成期で、電子の波動性という概念が導入され、経験的に見出されていた化学原理に対して理論的な裏付けがなされた。同時に新しい化学原理も多数発見された。しかし、大胆な近似のため本質を表せる化学現象が限られていた。
第 2 期は理論化学の確立期ともいうべき時期で、量子化学理論の基盤が確立され、コンピュータを用いた計算化学も産声を上げた。高精度な励起状態計算である SAC-CI 法(1978 年)や今日 gold
standard と評される CCSD(T)法(1989 年)もすでにこの時期に開発されている。第 3 期はまさに発展期で、コンピュータの飛躍的な発達も伴って、実用性の高い量子化学計算法の開発が行われた。
この時期のトレンドは、大規模化・高精度化・長(実)時間化であり、大ブレークしたのが密度汎関数理論(DFT)である。本講演者の研究経歴のほとんどがこの時期に属し、これらのトレンドに沿った研究も数多く行ってきた(表 1)。本講演者は、化学的理解には解析も重要と考え、解析手法の発展も行ってきた(表 1)。
これらの発展により理論化学が数値的に予言できる領域が広がったのは疑う余地もない。実際、分子構造・熱力学量・分光学的定数・物性値は、時には実験以上の精度で計算できる。しかし、高度化した計算結果をもとに、普遍的な化学原理を導くことが困難になったのも事実である。本講演では、まず化学原理について、具体的な研究例を紹介し、考えてみたい。表 2 は、本講演者らが行った応用研究のうち、化学原理まで踏み込んで議論できたものを示している。応用研究(1)~(3)はいずれも実験研究者から話題提供を受けて、それに応える形で進めた研究である。(1)では、従来の超共役とは異なる励起状態に特有なπ*-σ*超共役を見出すことができ、実験研究者に
表 1．講演者らの開発した理論・方法論
方向性
大規模化

理論・方法論
分割統治法

高精度化

核・電子軌道
理論

高精度化

内殻励起
汎関数
局所応答
分散力法
分子軌道
内挿法
実時間発展法
エネルギー
密度解析

高精度化
長時間化
実時間化
解析手法

略称
DC-HF, DC-DFT, DC-DFTB,
DC-MP2, DC-CCSD, DC-CCSD(T),DC-SACCI, DC-TDDFT, DC-GF
NOMO/HF, NOMO/MP, NOMO/CI,
NOMO/CC, TF-NOMO, TRF-NOMO
ECG-NOMO/HF, ECG-NOMO/MP2
CV-B3LYP, CVR-B3LYP, OS
LRD
LIMO, LSMO
RT-TDHF, RT-TDDFTEDA, BEDA, NAO-EDA, Gird-EDA,
IT-EDA, Interaction-EDA, ETS

概要
線形スケーリング法

総説
3,4

非断熱分子理論

5,6

高精度な内殻励起状態計算
法
第一原理的な分散力計算法

7

AIMD シミュレーションにおける SCF 収束性改良法
電子状態の実時間発展法
エネルギー分割による解析8表 2．講演者らの行った応用研究と関連する化学原理
研究テーマ(1) 励起状態におけるメチル基内部回転(2) ジベンゾスベレン
化合物の緩和過程(3) 表面吸着した環状
水 6 量体の分子振動(4) 高対称性分子の
励起状態

実験研究
超音速ジェット電子スペクトル分光法（藤井正明氏）
ピコ秒・ナノ秒過渡吸収スペクトル（高橋博彰氏）
赤外反射分光（伊藤正晴氏）
可視・紫外吸収スペクトル

化学原理
Hammet 則
π*-σ*超共役
スピン－軌道相互作用
El-Sayed 則
振電相互作用IRAS 活性非垂直振動
縮重系励起の対称則

総説

9,10

も受け入れられている。一方、応用研究(4)の「縮重系励起の対称則」は、数多くの励起状態計算を行うなかで、本講演者らが見出した化学原理である。理論化学の強みは、対象を選ばないことである。「縮重系励起の対称則」は、CO2, N2 などの小分子からアセチレン
やベンゼンなどの有機分子、そして金属錯体、金属ク
ラスター、フラーレンなどでも成立する。是非、高対
称性分子の励起状態を扱う実験・理論研究者には、この化学原理を理解し、利用してもらいたい。
理論化学が直面するもう一つの問題は、周期表の全元素を高精度に取り扱えるかということである。これには、相対論効果と電子相関効果を同時に考慮する必要がある。解くべき基礎方程式も、非相対論的な 図 1. IODKH/IODKH ハミルトニアン（LUT
Schrödinger 方程式から相対論的な Dirac 方程式にシフ スキーム有無）に対する従来法および DC 法
トしなければならない。実用的な手法には、もちろん による HF, MP2, CCSD 計算の計算時間。計算
大規模化も欠かせない。本講演者らは、最近、この問 系は鎖状(HF)n 分子(n = 2, 3, …, 50)。計算機
は、Hexa Core Xeon/3.33 GHz プロセッサーの題にも取り組み、一定の方向性が見えてきた(図 1)。
今後の 20 年が理論化学にとってどのような時期になるか、充実期なのか、変革期なのか、現時点では不明である。しかし、まだまだ理論化学が解決すべき問題は山積している。化学原理の探究と相対論的量子化学へのパラダイムシフトというのは全く別のスタンスに見えるが、両者を同時に進めるメリットはあると信じている。
[1] 「巨大分子系の量子化学法」、中井浩巳、CSJ カレントレビュー・シリーズ第 8 巻『巨大分子系の計算化学－
超大型計算機時代の理論化学の新展開』(日本化学会編)、pp. 52-60、(化学同人、2012 年).
[2] 「理論化学における理論の革新」、中井浩巳、別冊化学『化学のブレークスルー【理論化学編】』(日本化学会編)、pp. 66-72、(化学同人、2010 年).
[3] “Divide-and-conquer approaches to quantum chemistry: Theory and implementation”, M. Kobayashi, H. Nakai, pp.
97–127 in ‘Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications’, R. Zalesny, M.
G. Papadopoulos, P. Mezey, J. Leszczynski (Eds.) (Springer, 2011).
[4] 「分割統治(DC)電子状態計算プログラムの GAMESS への実装」、M. Kobayashi, T. Akama, H. Nakai, J. Comput.
Chem. Jpn., 8, 1-12 (2009).
[5] “Nuclear orbital plus molecular orbital (NOMO) theory: Simultaneous determination of nuclear and electronic wave
functions without Born-Oppenheimer approximation”, H. Nakai, Int. J. Quant. Chem., 107, 2849-2869 (2007).
[6] 「原子核と電子の波動関数を同時に決定するための non-Born-Oppenheimer 理論：NOMO 理論」、中井浩巳、星野稔、Mol. Sci. (分子科学), 1, A0010 1–22 (2007).
[7] “Description of core ionized and excited states by density functional theory and time-dependent density functional
theory”, Y. Imamura, H. Nakai, pp. 275–308 in ‘Quantum Systems in Chemistry and Physics: Progress in Methods andApplications’ Progress in Theoretical Chemistry and Physics’, B 26, K. Nishikawa, J. Maruani, E. J. Brändas, G.
Delgado-Barrio, P. Piecuch, (Eds.) (Springer, 2012).
[8] “Energy Density Analysis of Kohn-Sham DFT Method and Its Applications”, T. Baba, Y. Yamauchi, Y. Kikuchi, Y.
Kurabayashi, H. Nakai, Bull. Soc. Discrete Variational Xα, 18, 7-19 (2005).
[9] “Rules for excited states of degenerate systems: Interpretation by frozen orbital analysis”, H. Nakai, pp. 363–395 in‘Advances in the Theory of Atomic and Molecular Systems: Conceptual and Computational Advances in Quantum Chemistry’
Progress in Theoretical Chemistry and Physics, P. Piecuch, J. Maruani, G. Delgado-Barrio, S. Wilson (Eds.) (Springer, 2009).
[10] 「化学原理の発見：縮重系励起の対称則」、H. Nakai, J. Comput. Chem. Jpn., 11, 1-12 (2012).