4P104
周期境界条件の ONIOM 法への適用（上智大院・理工）○小林 理，南部 伸孝

Application of Ewald Sum for ONIOM method(Graduate school of Sophia Univ.) ○Osamu Kobayashi, Shinkoh Nanbu【序】諸熊らにより考案された ONIOM (our Own N-layer integrated molecular Orbitals andmolecular Mechanics)法は量子化学計算と分子力場法を統合し、巨大系に対する計算を可能にする手法として知られている[1]。ONIOM 法は系の全エネルギーを以下のように分割する。

EONIOM = EMM,real + EQM,model − E MM,model
ここで EMM,real 、 EQM,model 、 EMM, model はそれぞれ系全体に対する molecular mechasnics (MM)法によるエネルギー、model 系に対する quantum mechanics (QM)法によるエネルギー、model 系に対する MM 法によるエネルギーである。この方法は mechanical embedding 法とよばれる。また同じく諸熊らにより QM 層と MM 層の相互作用を適切に取り入れるために、MM 層の電荷を点電荷として QM 層の計算に取り込む electronic embedding 法が開発された。この手法では MM層の原子と QM 層の電子の相互作用は、1 電子演算子で表現される。
MM が十分に大きな系であるとき、 EMM, model の計算を有限の大きさで近似して行うと、表面が生じることに起因し問題が生じる。このような問題は例えば溶液や結晶においてみられる。これは周期境界条件(PBC)の導入により解決される。Coulomb エネルギーは原子間距離に対して収束が非常に遅いことが知られており、PBC においては 2 つの収束の速い級数に分割し合計する
Ewald sum 法や計算を効率化した particle mesh Ewald sum 法(PME)が広く用いられている。
我々はこれまで非断熱現象を取り扱うために ONIOM 法を Tully’s surface hopping (TSH)法と組み合わせた QM/MM-TSH 法による分子動力学(MD)計算を行ってきた[2]。しかし、従来の方法では有限系の取扱いに限られており、系の周囲の環境が適切に取り込まれていない可能性があった。そこで
我々は ONIOM 法に Ewald 法を組み合わせた ONIOM-Ewald 法を開NH2非断熱現象発し、さらにそれを QM/MM-TSH 法に適用した。本発表では非断熱現象により励起状態における cis 型から trans 型への異性化が知られている 2,4-pentadiene-iminium cation (protonated Schiff base 3,PSB3)のメタノール中の MD 計算を行い、本手法の効果を調べる。
【方法】シミュレーションを行う系の構造と分子力場は以下の手順に

NH2

図 1. PSB3 における非断熱的な構造変化より作成した。はじめに PSB3 を Gaussian 09 を用いて構造最適化を行った。構造最適化は密度汎関数法により行い、密度汎関数は B3LYP、基底関数は cc-pVDZ を用いた。引き続き AMBER9プログラムパッケージに付属する ANTECHAMBER を用い restrained electrostatic potential(RESP)電荷を決定した後、メタノールを PSB3 の周囲に配置した。得られた系の分子力場はgeneral AMBER force field (GAFF)を用いた。作成した系は AMBER9 の Sander モジュールを

使い 1 atm, 300K での NPT アンサンブルで 1 ns 時間発展させ、次の計算に用いた。
構造変化を調べるために、QM/MM-TSH 法を用い励起状態 MD シミュレーションを行った。
シミュレーションの開始時点において、各原子の初期速度はボルツマン分布に従いランダムに発生させ、PSB3 の電子状態は S1 にあると仮定した。シミュレーションは NVE アンサンブルで行い、タイムステップは 0.50 fs とした。遷移確率の計算には朱-中村理論を用いた。各時刻においてエネルギーと各原子にはたらく力は Gaussian 09 を用い ONIOM-Ewald 法に従い求めた。
PSB3 を QM 層に、その他の分子を MM 層に割り当てた。
8ONIOM 法は単位セルに関して electronic embedding 法を6用い、外側の像に関しては mechanical embedding 法を用いた。QM 層の計算は Molpro 2012 を、MM 層の計算は当研究室で開発したオリジナルプログラムを用いた。MM 層の計算で用いた分子力場は先に初期構造を得るために用Potential Energy / eV

この割り当てはシミュレーション中変更していない。

(a)

4S120

いた分子力場と同じものを用いた。また QM 層の計算は-2CAS(6e,6o)SCF 法を用い、基底関数には MIDI-bang を用-4-4S000いた。比較のために周期境界条件を適用しない electronic

200200400400

embedding 法による MD 計算もあわせて行った。この計8ONIOM-Ewald 法を用いたときと同様である。いずれの場6対して行い、S0 状態への非断熱遷移が確認された後はシミュレーションを停止した。
【結果と考察】はじめに分子の各層への割り当てを固定したことの妥当性は、このシミュレーションにおいては単位セルから外部への溶媒分子の拡散はほとんど観察されな

Potential Energy / eV

算の詳細は外場の取り扱いを変更した以外は合もシミュレーションは現在 100 本のトラジェクトリに60060080080010001000time / fs(b)4S12S00-2-40

かったことより支持される。これはシミュレーション時間

200400600

800

1000time / fsが溶媒の拡散速度に対して十分に小さかったためである。
図 2 にシミュレーションに用いた全てのトラジェクトリに図 2. PSB3 の S0 状態及び S1

ついて PSB3 のポテンシャルエネルギーの時間変化を示し状態のポテンシャルエネルギーた。図 2 において(a)は 1000 fs までグラフが続き非断熱遷(a) ONIOM-Ewald 法による移せずに S1 状態を保つトラジェクトリがあることを示す(b) cut-off による一方、(b)は 400 fs 付近においてグラフが途切れ全てのト(実線が実際の電子状態を表す)ラジェクトリにおいてこの時刻までに S0 状態に非断熱遷移したことを表している。これは ONIOM-Ewald 法を用いた場合、非断熱遷移が遅くなることを意味している。Muñoz-Losa[3]らは、非断熱遷移が無輻射遷移に要する時間は溶媒の運動に影響されると考察しており、我々の結果は Muñoz-Losa の結果を支持している。
[1] T. Vreven et. al. J. Chem. Theory Comput., 2006, 2, 815-26[2] W. C. Chung et. al. J. Phys. Chem. B, 2012, 116, 8009–8023[3] A. Muñoz-Losa et. al. J. Chem. Theory Comput., 2011, 7, 4050-4059