3E14

分子数適応階層型 QM/MM-MD 法による溶質-溶媒間相互作用の改善：グリシン水溶液への適用(名大院・情報科学 1, JST-CREST2) 竹中規雄 1,2, 北村勇吉 1, 小谷野哲之 1,2, 長岡正隆 1,2

An Improvement in Solute-Solvent Interactions via theNumber-Adaptive Multiscale QM/MM-MD Method: Application
to Glycine Aqueous Solution(Graduate School of Information Science, Nagoya University1, JST-CREST2)
Norio Takenaka1,2, Yukichi Kitamura1, Yoshiyuki Koyano1,2, Masataka Nagaoka1,2【序】量子力学的(QM)/分子力学的(MM)分子動力学(QM/MM-MD)法は、凝集分子系において、着目する溶質分子の電子状態変化を考慮しながら、多数の溶媒分子群の運動も同時に取り扱うことが可能であり、計算精度と計算時間のバランスが良いため非常に広く用いられている。しかしながら、その計算精度は、溶質の QM 計算の精度だけではなく、溶質-溶媒(QM-MM)間相互作用の精度にも強く依存しており、主に 2 つの問題点が良く知られている。一つ目は、QM-MM 間の静電相互作用が近距離で不正確になる問題であり、これは QM 溶質の電子が MM 電荷へと引き寄せられ、その電子配置が不正確になることに起因する[1]。二つ目は、Lennard-Jones 関数等を用いて経験的に取り扱われる QM-MM 間の非静電相互作用の記述の問題であり、一般的にはGAFF 力場等の汎用的な力場が良く用いられているものの、QM/MM-MD 法においても適切であるとは限らない。このような問題点を解消するため、溶質-溶媒(QM-MM)間の相互作用を“経験的”に改善する方法がこれまでに幾つか提案されてきた[1-3]。しかしながら、これらの方法では、適切な参照値を必要とするため汎用性に乏しい点と化学反応過程における中間状態での取り扱いが困難である点が大きな欠点である。そこで、これらの問題点を“非経験的”に解消するため、本研究では、溶質近傍の溶媒もQM 領域へ含めた計算が実行可能な分子数適応階層型 QM/MM-MD 法[4]を新しく提案し、本手法をグリシン水溶液系へ適用してその有効性を検証した[5]。
【理論と方法】溶質近傍の溶媒も QM 領域へ含めた場合には、溶媒が QM 領域と MM領域に分割されるため、それらの間で溶媒分子の入れ替わりが起こり、エネルギー勾配が不連続になってしまう点が大きな問題である。このため、2 つの領域の間に遷移(T)領域を設けて、入れ替わる溶媒分子の勾配を連続的に変化させる必要がある。この際に、T 領域は、従来は溶質からの“距離”を用いて定義されていたが[6]、我々は計算上の優位性から、溶媒の“分子数”を用いた取り扱いを新しく提案した[4,5]。本手法
では、MD の各ステップで T 領域に含まれる NT 個の溶媒分子を 1 個ずつ QM 領域へ加えることにより構築される合計 NT+1 個の QM/MM ポテンシャル Va (a = 1,…,NT+1)を用いて、運動方程式とポテンシャルエネルギーVad が以下のように与えられる。
MNT 1d 2r
 a  r  t   V a  r  t  dt 2a 1

 NT 1 a    r  t    1
 a1
NT 1NT 1V ad  r  t      a  r  t  V a  r  t       a  r  t  V a  r  t  dra1r

a1(1)(2)

4

4SD (OC-Ow )

3.5
32

1.5

11

0.5

0.5

CαCc21.50OC2NAM (N-Hw )2.5g(r)g(r)SD (N-Hw )NAM (N-Ow )3NAM (OC-Hw )2.5SD (N-Ow )3.5SD (OC-Hw )NAM (OC-Ow )

HαN

O1C
H 3N

H1,2N01234567r (Å)(a) -COO 基(OC)-水分子(OW,HW)間1234567r (Å)(b) -NH3+基(N)-水分子(OW,HW)間

図 1 双性イオン型グリシン分子の周りの水分子の動径分布関数(SD: 標準的な QM/MM-MD 法、NAM: 分子数適応階層型 QM/MM-MD 法)

ここで、M と r は系の質量と座標のベクトルであり、σa は分割 a の重み関数である。
本手法は、溶質のみを QM とする標準的な QM/MM-MD 法に比べると計算コストは高いものの、式(1)における NT+1 個の勾配は各ステップで並列に独立して計算が可能であるため、実質的には最大の QM 領域(a = NT+1)を持つ QM/MM 計算のコストと同程度まで計算コストを削減することが可能であり、これは本手法の大きな利点である。
【結果と考察】本研究ではテスト系として、水溶液中の双性イオン型グリシン(NH3+CH2COO-)へ本手法を適用し、水和構造の解析と Natural Population Analysis による電荷分布の解析を行った[5]。本解析では、QM 計算の理論レベルに B3LYP/6-31G(d)法、溶媒水分子には TIP3P モデルを採用した。この際に、溶質のみを QM とした場合(標準的な QM/MM-MD 法)と溶質の周りの第 1 水和圏の水分子も QM 領域へ含めた場合(分子数適階層型 QM/MM-MD 法)の 2 つの計算結果を比較した。図 1 にグリシンの(a) -COO-基および(b) -NH3+基の周りの水分子の動径分布関数をそれぞれ示す。標準的な QM/MM-MD 法では、全系を QM とする Car-Parrinello MD (CPMD)法による結果[7]に比べて、動径分布関数の第 1 ピークは両方とも高くなり、結果としてこれらの水和数は実験値よりも過大に見積もられた。また、水溶液中での電荷分布の解析によると、
QM 溶質の電子と MM 電荷間でパウリ反発が欠如しているため[1]、-COO-基と-NH3+基は過剰に分極しており、これが溶質-溶媒間相互作用を過大評価する原因であることが判った。一方、分子数適応階層型 QM/MM-MD 法は、CPMD 法による水和構造や水和数の実験値を精度良く再現した。したがって、溶質近傍の溶媒分子を QM 領域へ含めることにより、溶質-溶媒間相互作用を十分に改善することが可能であり、今後、凝集分子系の化学反応の自由エネルギー解析[8]等への適用が大いに期待できる。
[1] A. Laio, J. VandeVondele, U. Rothlisverger, J. Chem. Phys., 116, 6941 (2002).
[2] Y. Tu, A. Laaksonen, J. Chem. Phys., 113, 11264 (2000).
[3] (a) N. Takenaka, Y. Koyano, M. Nagaoka, Chem. Phys. Lett., 485, 119 (2010). (b) Y. Koyano, N.
Takenaka, Y. Nakagawa, M. Nagaoka, Bull. Chem. Soc. Jpn., 83, 486 (2010).
[4] N. Takenaka, Y. Kitamura, Y. Koyano, M. Nagaoka, Chem. Phys. Lett., 524, 56 (2012).
[5] N. Takenaka, Y. Kitamura, Y. Koyano, M. Nagaoka, J. Chem. Phys., 137, 024501 (2012).
[6] (a) T. Kerdcharoen, K. R. Liedl, B. M. Rode, Chem. Phys., 211, 313 (1996). (b) T. Kerdcharoen, K.
Morokuma, Chem. Phys. Lett., 355, 257 (2002). (c) A. Heyden, H. Lin, D. G. Truhlar, J. Phys. Chem. B,111, 2231 (2007). (d) R. E. Bulo, B. Ensing, J. Sikkema, L. Visscher, J. Chem. Theor. Comp., 5, 2212
(2009).
[7] J .Sun, D. Bousquet, H. Forbert, D. Marx, J. Chem. Phys., 133, 114508 (2010).
[8] N. Takenaka, Y. Kitamura, Y. Koyano, T. Asada, M. Nagaoka, Theor. Chem. Acc., 130, 215 (2011).