1E16

構造最適化 Elongation 法の開発と局所振動解析（九大院・総理工1, SCNU2, JST-CREST3）○青木百合子1,3、Liu Kai1, Yan Yun-an1, 折本 裕一1, Feng Long Gu2,3
Development of elongation optimization method and local vibrational analysis(Kyushu Univ.1, SCNU2, JST-CREST3)

Yuriko Aoki1,3, Liu Kai1, Yan Yun-an1, Yuuichi Orimoto1, Feng Long Gu2,3【序】量子化学計算は飛躍的に発展を遂げたものの、未だ生体高分子や固体などの巨大系に対しての高精度かつ高速計算は依然として困難である。現在までに我々は、局在化と相互作用を繰り返しながら高分子の重合反応を計算機上で実現することにより、効率的に電子状態を合成するための新規な方法を開発してきた。近年、一次元系を仮定したこのElongation(ELG)法を、超高精度で二次元・三次元系に適用可能となるように発展させ、高速かつ厳密にバルク系の電子状態をも計算できる 3D-Elongation(3D-ELG)法に展開した。大規模系における化学反応を扱う際に、例
え化学反応部位が局所的現象であっても、従来法によると全系の構造最適化や振動解析が必要であり、系の増大とともに扱いが困難となってくる。そこで今回、本方法に効率的局所構造最適化法(ELG-OPT法)を組み込んだ結果、構造最適化等が高効率で可能になっただけでなく、従来法では到達し得ない安定構造に到達できる可能性が示唆され、さらに局所振動解析法(ELG-VIB法）を組み込むことにより効率的な反応部位の振動解析が可能となった。例えば非線形光学特性を算出するには、全エネルギーの電場に対する高次微分が必要になり、例え 10-6 a.u.程度の小さな誤差でも、特にγ値においては 10～20％の誤差を生む。同様な状況が、全エネルギーの座標に対する二次微分を要する振動解析にも当てはまり、僅かなエネルギー誤差も振動数への大きな誤差へとつながる。よって特に遷移状態探索を局所的にかつ全系を扱ったときと同じ精度で実行するには、全系の性質を取り込んだ上での反応部分の正確な振動構造を抽出するノウハウが必要となる。本研究では、高分子反応における触媒活性の解析に用いることを目指し、まずElongation法による局所構造最適化法の組込を行ない、さらに局所振動解析が可能となるように発展させるとともに、効率的な局所遷移状態探索に拡張することを目的とする。
【方法】本方法では、系の伸長の過程で、反応スポットにおいて攻撃分子(M)との相互作用が必要な領域をActive領域とし、予めそこに局在化しておいた領域局在化軌道(RLMO)のみがMと相互作用し、それ以外の領域をFrozen領域として、そこに局在化している凍結領域局在化軌道(FrozenRLMO) を保存していく。さらに、Active RLMO基底のFock行列に影響を及ぼさないAO基底Fock行列について、寄与しないことが分かっている二電子積分の計算を最初から回避し、反応部の影響を受けない構造パラメータを固定することによって勾配の計算を高速化した（Cutoff-ELG-OPT法）。
ここで、全系計算におけるエネルギー勾配の各項において、電子密度はELG法で得られたものは従来法と同じであるのでそのまま利用できるが、重なり積分の勾配の係数に現れる軌道エネルギーの項は、ELG法では非対角行列となることが従来のエネルギー勾配との異なる点である。
一方、遷移状態探索は化学反応解析において重要であるが、大規模系の場合は全系に対する振動解析は現実的ではない。しかし反応部分のみの力の定数(Hessian)行列が正確に分かれば反応に関わる遷移状態探索は可能である。ただし反応中心の部分 Hessian 内で遷移状態を探索したのでは、そこだけくり抜いて計算したクラスター計算と同じである。ELG 法では、反応に関与しない部分の Hessian を計算することなく、その部分の影響も含んでいる Active 領域（B＋M）に対するHessian を計算するため、Frozen と Active の境界領域の Hessian も正しく表すことができる（ELG-VIB 法）。本方法では、系全体の電子密度行列を、Active 領域と Frozen 領域からの寄与に分類できるので、Active 領域のみの有効 Hamiltonian から部分 Hessian 行列が得られるが、この際、Active 領域と Frozen 領域間の相互作用は多重極展開によって見積もる。

1(nuc)(BM)  (BM)(BM)(nuc)(e-nuc)E ( eff ) = VBM+ ∑ D(BM)µν D ρσ ( µν ρσ ) +VA,BM + VA,BM,BM + ∑ D µν H µν
2 µνρσµν (BM) = 1 ∑ D(A) ( µν ρσ ) +H (BM)Hµνρσµν2 ρσ
これにより、反応に関与しない部分の Hessian を計算することなく、その影響下で反応部分のHessian を正しく計算し、小行列の対角化を経て効率よく局所振動モードを得ることができる。
【結果】Elongation-構造最適化(ELG-OPT)法を適用したところ、ほぼ正確に従来の構造最適化(CONV-OPT)法を再現するか、あるいはより安定な構造を見出すケースがあることが分かった。
ELG-OPT 法によって得られた最適構造がCONV-OPT法による結果と一致する場合は、全エネルギー誤差は～10-8a.u./atomを達成しているが、図１のDNAの例のように、ELG-OPT法によって到達し得た最適構造の方が安定な場合がある(～－8×10-4 a.u./atom安定)。
その構造でのCONV法（図 1 表中「CONV single_point」）とは完全に一致しているため（～10-8a.u./atom)、マイナスの差はELG法から派生したものではなく、従来法が準安定構造で最適化を停止したことからくる。これは、同じ初期構造から最適化を行なった場合に、従Red: CONV : Blue: ELG

RMSD: 0.25来法に比べてELG-OPT法ではよりグローバルな安定構造探索が可能

Total energy by
ELG-OPT (a.u.)

Diff. (a.u./atom)であることを示唆しており、従来の構造最適化では到達し得ない安-27293.8266632312

定構造を見出すケースと考えられる。その理由としては、全系をまCONV single_point-27293.8266884810

ともに最適化する従来法では初期構造に強く依存し、一旦Local逐次的最適化を行うため、全系の最適化が容易になることが考えられる。
次に、ELG-VIB法によって得られた振動モードの

直線状 水分子鎖
H2O…H(H2O…H2O…H2O…H2O…H2O…H2O…H2O…H2O…H2O…H2O…H2O2O)20
A化されていることは確認済みである。いずれも反応末端部に起因したモードであるが、水分子鎖でH2O…H2O

B

M

ELG (cm-1)

4375.06

4377.994380.314385.824388.42

CONV (cm-1)

4387.104387.33

4390.314391.184398.05

直線状 ポリグリシン (-CO-N-C-)20

比較を３つの系に対して図２に示す。局所Hessian
作成に要する演算時間は従来法に比べて相当高速

3.85E-08

図１ DNA の ELG-OPT およびCONV-OPT による構造最適化

minimum に 落 ち 込 む 確 率 が 高 い の に 対 し て 、
ELG-OPT法では、各部分空間で形状記憶をしながら-8.10E-04

ELG-OPT法が安定

A

M

B

ELG (cm-1)

4054.984055.274055.624056.324057.51

CONV (cm-1)

4058.98

4060.244058.214058.854060.86

ポリアラニン、Random β-sheet (N=16)は10cm-1、化学結合系では数cm-1以内で一致している。しかし、これらは最適構造での比較ではない。

AB

MM

M

ELG (cm-1)

4244.26

4238.764211.434115.474096.07

上記最適化法と組み合わせることにより、エネル

CONV (cm-1)

4243.534237.18

4210.364114.404094.73

ギー極値での振動モード解析を行う予定である。

図２ 水分子鎖、ポリグリシン、ポリアラニンのELG-VIB 法による振動モードの比較