2E14交換対称性を利用した QM/MM 法の境界条件（原子力機構 1、サッサリ大 2）○志賀 基之 1、マシア マルコ 2

QM/MM (quantum mechanics/molecular mechanics) 法は、着目する部分系(QM 領域)を第一原理的な電子状態計算で扱い、その周囲の部分系(MM 領域)を経験的な古典力場計算で扱うという、ハイブリッド法である[1]。全電子計算と比べて計算量が格段に減るので、溶液や蛋白質など複雑多体系を扱うのに適している[2-5]。ただし、境界条件の決め方など、計算法としては未熟で発展途上の面もある。特に最近では、分子動力学法など、分子シミュレーションとの組み合わせにおいて、計算の信頼性や効率化が議論の対象になっている[6-10]。
本研究で取り上げる問題は、QM 領域と MM 領域の間を行き交う分子が存在するときの不都合である。例として、希薄なイオン水溶液系を考えよう。いま、イオンと第一水和圏の水分子に間に働く複雑な相互作用を第一原理計算で扱いたいので、イオンとイオンに最も近いn 個の水分子を QM 領域、残りの m 個の水分子を MM 領域に設定し、m>>n とする。この条
件で分子動力学シミュレーションを実行すると、水和圏内の水とバルクの水が自然に入れ替わってしまう。もしこれを放置すれば、やがてイオンの周りには QM の水分子はいなくなり、MM の水分子により取り囲まれるだろう。m>>n の場合、交換したほうがエントロピー的に有利なので、イオンと水の相互作用の強さとは無関係に、時間が経てば交換はいつか必ず起こる。しかし、これでは第一水和圏をきちんと扱おうとした目的にそぐわない。
最も安易な解決法は、水和圏内の水とバルクの水の入れ替わる瞬間に、QM と MM のラベルを入れ替えることである。しかし、ふつう QM/MM ポテンシャルは QM 原子と MM 原子の座標交換に対して対称性が崩れているので、そこでポテンシャルが不連続になり、分子動力学シミュレーションに支障をきたす。これを防ぐのに、ポテンシャルの段差を和らげる方法がいくつか提案されている[11-16]。
本研究では、もう一つの解決法として、外部から荷重ポテンシャルを加えて交換を抑制する方法を提案する[17]。荷重ポテンシャルは人工的なものだから、いかなる場合も動的な情報は損なわれる。しかし、時間に依存しない熱力学量（イオン水溶液では、動径分布関数、水和エネルギーなど）を全く変えないようにすることが可能である。これは、厳密な定式化によれば、QM 原子と MM 原子の座標交換で得られるあらゆる配置に対して存在確率を適切に決め、それに応じて荷重ポテンシャルを与えればよい[17]。しかし、実際の QM/MM 分子動力学計算では、次のように最も簡略化された“単交換近似”でも、かなりうまくいくことがわかった。これには、まず、原点（イオン水溶液ではイオンの位置とする）から最も遠い QM
原子 i への距離 ri と、原点から最も近い MM 原子 j への距離 r j の差 rij  rj  ri を計算する（図１）。次に、 rij を変数とした Fermi 関数の対数を荷重ポテンシャルとしてV bias  k BT log f ij ,f ij 

11  exp rij と与える。ここで、k B は Boltzmann 定数、T は温度であり、 f ij と f ji はそれぞれ配置(i,j)と(j,i)の現れる確率因子を表していて、和が 1 である。最後に、 V bias と、その力   R V bias をQM/MM 計算において加えるだけである。これは実装が簡単で、計算負荷も極めて少ない。

図１：原子の座標交換。●は QM 原子、○は MM 原子を表す。
ここで、指数  は QM 原子と MM 原子の分離しやすさを制御する入力パラメータである。

 を無限大にすれば rij が常に負となり、完全に分離される。これは都合がよい反面、荷重ポテンシャルが急峻になってしまう問題がある[18]。したがって、実際の分子動力学計算では

 を大きくしすぎないようにすべきである。
以上の方法は、荷重ポテンシャルの原点の定義を変えれば、界面など、他の境界にも応用できる。また、QM/MM 法に限らず、ONIOM 法などハイブリッド法全般に使えるものである。
本発表では、多重時間ステップ法[6,7]と組み合わせた効率のよい QM/MM 分子動力学計算の具体例を紹介し、その有用性を示す予定である。
[1] A. Warshel, M. Levitt, J. Mol. Biol. 103, 227 (1976). [2] H. Lin, D. G. Truhlar, Theor. Chem. Acc.
117, 185 (2007). [3] H. M. Senn, W. Thiel, Angew. Chemie Intl. Ed. 48, 1198 (2009). [4] H. Hu, W.
Yang, Ann. Rev. Phys. Chem. 59 (2008). [5] M. Boeckmann, D. Marx, C. Peter, L. D. Site, K. Kremer,N. L. Doltsnis, Phys. Chem. Chem. Phys. 13, 7604 (2011). [6] M. Shiga, M. Tachikawa, Mol. Simul.
33, 171 (2007). [7] T. K. Woo, P. Margl, P. R. Blochl, T. Ziegler, J. Phys. Chem. A 106, 1173 (2002).
[8] H. Takahashi, Y. Kawashima, T. Nitta, N. Matsubayasi, J. Chem. Phys. 123, 124504 (2005). [9] T.
Yamamoto, J. Chem. Phys. 129, 244104 (2008). [10] T. Janowski, K. Wolinski, P. Pulay, Chem. Phys.
Lett. 530, 1 (2012). [11] T. Kerdcharoen, K. Morokuma, Chem. Phys. Lett. 355, 257 (2002). [12] B. M.
Rode, T. S. Hofer, B. R. Randolf, C. F. Schwenk, D. Xenides, V. Vchirawongkwin, Theor. Chem. Acc.
115, 77 (2006). [13] A. Heyden, H. Lin, D. G. Truhlar, J. Phys. Chem. B 111, 2231 (2007). [14] R. E.
Bulo, B. Ensing, J. Sikkema, L. Visscher, J. Chem. Theory Comput. 5, 2212 (2009). [15] N. Takenaka,Y. Kitamura, Y. Koyano, M. Nagaoka, J. Chem. Phys. 137, 024501 (2012). [16] N. Bernstein et al.
Phys. Chem. Chem. Phys. 14, 646 (2012). [18] M. Shiga and M. Masia, J. Chem. Phys. 139, 144120(2013). [17] C. N. Rowley and B. Roux, J. Chem. Theor. Comput. 8, 3526 (2012).