3E10
ONIOM 法と Microiteration 法を用いた反応経路自動探索法の巨大系への応用: イソペニシリン N 合成酵素への応用（北大院理 1，京大福井謙一研究セ 2，豊田理研 3，エモリー大化学 4）前田 理 1，阿部 恵里
花 2，畑中 美穂 2，武次 徹也 1，大野 公一 3，諸熊 奎治 2, 4

Exploring Potential Energy Surfaces of Large Systems with Automated Reaction
Path Search Methods in Combination with ONIOM and Geometrical
Microiteration: A Case Study on Isopenicillin N Synthase(Hokkaido University1, Kyoto University2, Toyota Physical and Chemical Research Institute3, Emory
University4) Satoshi Maeda1, Erika Abe2, Miho Hatanaka2, Tetsuya Taketsugu1, Koichi Ohno3, andKeiji Morokuma2,4【序論】化学反応は、反応経路を量子化学計算によって求めることによって理論的に解析できる。反応経路の計算において最も重要なのが、遷移状態(TS)の構造を求めることである。
TS 構造が求まれば、そこを通過する反応経路を自動計算できるからである。通常、TS は構
造最適化計算によって求める。近年、量子化学計算手法と計算機技術の発展により、数十原子からなる系の TS 構造はルーチン的に求めることができる。一方、ONIOM 法などの QM/MM
ハイブリッド計算法の開発により、反応中心を取り囲む多数の溶媒分子や巨大なタンパク質を考慮したエネルギー計算が可能になった[1]。そこで、巨大系の構造最適化計算のために
Microiteration 法が導入され[2]、数千または数万原子からなる酵素の中のローカルな反応中心
における化学反応について、全(数万)自由度を考慮した TS 構造最適化が可能となった[3]。
構造最適化で得られる構造は初期推定に依存する。これまで、初期推定は計算者の経験と直観に基づいて行われてきたが、複雑な反応機構について全ての TS 構造を推定することは難しい。また、反応機構が全く未知である場合には、関係する TS 構造を見落としてしまう危険もある。そこで我々は、経験や直観に頼らず系統的に反応経路を自動探索できる方法を開発してきた。超球面探索法では、安定構造から出発する異性化・分解経路を自動探索できる[4]。
人工力誘起反応法では、二つ以上の反応物の組から出発する生成経路を自動探索できる[5]。
これらを組み合わせれば、生成・異性化・分解、全てを扱える。本研究では、ONIOM 法および Microiteration 法と組み合わせ、イソペニシリン N 合成酵素(IPNS)への応用を行った。
【方法】Microiteration 法では、巨大分子を反応中心とそれ以外に分割し、反応中心原子を動かす度に、反応中心の構造を固定したエネルギー極小化(Microiteration)を実行する。これにより探索空間が大幅に減少し、効率的に TS 構造を探索できる。また、Microiteration を併用した探索で得られる TS 構造は、巨大分子の全自由度を考慮した TS 構造と完全に一致する。
本研究では、反応中心を B3LYP 法、それ以外を AMBER 力場で扱う QM/MM-ONIOM 法によってエネルギー勾配を計算した。IPNS には（タンパク内に存在する基質と水分子も含めて）合計 5368 個の原子が含まれ、その中の 65 原子を反応中心とした。エネルギー勾配計算および Microiteration は Gaussian09 プログラムで行った。反応中心の構造変化は、我々が開発してきた GRRM プログラム[6]に実装されている超球面探索法および人工力誘起反応法で扱った。

【結果】テストとして、水 100 分子に囲まれた HCHO 分子の異性化・分解経路の超球面探索法による自動探索[7]、および、水 300 分子に囲まれた HCHO + H2C=CH-OH 反応経路の人工力誘起反応法による自動探索[8]を行った。それぞれの計算を気相中で行った場合と比較し、気相中と同様の計算コストで数百原子系の計算を行うことができた。つまり、Microiteration 法が非常にうまく働いていることが分かった。
IPNS の反応中心には Fe 原子が存在し、Fe への基質の配位構造が X
線で観測されている。そこで、図 1 に示す構造からスタートし、超球

図 1、IPNS の反応中
心(7 重項)の構造。

面探索法を、障壁の低い経路のみを探索する[LADD = 3]オプ
シ ョ ン [6] を 使 用 して適用したところ、Fe に配位した O2 の配位構造が変化する反応経路を自動探索
することができた。
一方、第一ステップは 2 分子反応である
ため、人工力誘起反
応法が効率面から有効である。通常の人
工力誘起反応法の応用では、反応物同士

図 2、人工力誘起反応法による反応の第一ステップに関係する安定構造、
TS 構造、および、ポテンシャル交差領域内最小エネルギー構造(MSX) [8]。
各構造のスピン多重度は上付きでそれぞれの構造ラベルに示す。

を様々な配向から人工力で反応させる。一方 IPNS では、Fe 原子への基質の配位構造が X 線で観測されているため、これを利用することで計算量を軽減した。つまり、X 線で得られている反応物構造からスタートし、基質の官能基・置換基同士を人工力で反応させた。図 1 に、
基質の官能基・置換基の分割の様子を黒丸で囲んで示す。これら全てのペアに対して人工力
誘起反応法を適用することにより、第一ステップに関係する全反応経路を探索した。このとき、異なるスピン状態間のポテンシャル交差点も自動探索した。得られた構造群を図 2 に示す。この中で最も障壁が低い経路は 51 → 5TS1/3 → 53 → 5TS3/4 → 54 である。これは CH2 グループの H 原子を Fe に配位した O2 が引き抜く反応で、実験または過去の計算[3]から示唆されている反応と一致する。つまり、TS 構造や中間体の予想はもちろんのこと、生成物に関する情報すら一切用いずに、反応物の情報のみから実際に起こる反応を予測できた、と言える。
[1] Chung, L. W.; Hirao, H.; Li, X.; Morokuma, K. WIREs Comput. Mol. Sci. 2012, 2, 327. [2] Vreven, T.; Frisch, M. J.;
Kudin, K. N.; Schlegel, H. B.; Morokuma, K. Mol. Phys. 2006, 104, 701. [3] Lundberg, M.; Kawatsu, T.; Vreven, T.; Frisch,M.; Morokuma, K. J. Chem. Theory Comput. 2009, 5, 222. [4] Ohno, K.; Maeda, S. Phys. Scr. 2008, 78, 058122. [5] Maeda,
S.; Morokuma, K. J. Chem. Theory Comput. 2011, 7, 2335. [6] Maeda, S.; Osada, Y.; Morokuma, K.; Ohno, K. GRRM, adevelopmental version, 2012. [7] Maeda, S.; Ohno, K.; Morokuma, K. J. Chem. Theory Comput. 2009, 5, 2734. [8] Maeda,
S.; Abe, E.; Hatanaka, M.; Taketsugu, T.; Morokuma, K. (submitted).

3E11

SAC-CI 法による C60 fullerene の励起状態計算（1 分子科学研究所、2 計算科学研究センター、3JST CREST）
福田良一 1,2,3、江原正博 1,2,3

Electronic excitations of C60 fullerene calculated
by using SAC-CI method（1Institute for Molecular Science, 2Reserach Center for Computational Science,3
Japan Science and Technology Agency CREST）
Ryoichi Fukdua1,2,3, Masahiro Ehara 1,2,3
その科学的及び実用的な重要性から、C60 フラーレンは、量子化学計算の主要なターゲットである。その物性の特徴は、励起状態を含む電子状態に起因すると考えられるが、C60 の励起状態に対する高精度 ab initio 計算は、これまで報告されていない。フラーレン類の紫外―可視 (UV-vis)
スペクトルの帰属には、半経験的分子軌道法と TDDFT が用いられてきたが[1]、パラメータや汎関数依存性に対する疑念を払拭するためにも、多体波動関数法による高精度 ab initio 計算が求められている。しかし、分子サイズと非局在化した電子状態が、高精度 ab initio 計算を困難にしている。さらに、紫外―可視光領域には、非常に多くの状態が存在するため、スペクトルの帰属には、多くの励起状態を効率良く扱える理論が必要である。C60 に関しては、Ih 点群を十分に活用す
れば、こうした困難は、乗り越えられるかもしれない。しかし、そのような計算プログラムは、
一般には普及しておらず、また、さらに進んだ応用計算を考えると、高い対称性を利用することは現実的ではない。
最近我々は、direct SAC-CI 法[2]が、大きな共役分子系の計算に有効であることを見出した[3]。
direct SAC-CI 法は、Gaussian 09 B01 以降で公開されているプログラムである。今回は、その開発版で、主に並列化により計算効率を高めた direct SAC-CI プログラムを用い、C2h 点群の下で、C60の励起状態の計算を行い、UV-vis スペクトルの帰属を行った。求めたのは、Ih では 61 の励起状態
で、縮退を考慮しない C2h 点群で実際に計算したのは 230 状態である。
SAC-CI 法 は ク ラ ス タ ー 展 開 に 基 づ く 理 論 で あ り 、 基 底 状 態 を 指 数 関 数 型 の SAC(symmetry-adapted cluster)波動関数
 SAC  exp   cI S I†   0 ,(1)

Iにより記述し、励起状態は SAC 波動関数からの配置換相互作用(configuration interaction: CI) により記述する。
p
 SAC-CI  d Kp RK†  SAC .








(2)

K

ここで、R, S は１及び２電子励起演算子である。SAC-CI 法では、摂動論により、エネルギー的に重要な R, S のみを考慮し、重要でない励起演算子を、はじめから計算に含めないことで、計算量を大幅に減らしている。また、エネルギーベースの選択法であるから、非局在化した電子状態にも適用可能な方法である。また、SAC-CI 波動関数を求める永年方程式は、固有値問題であり、多数の状態を一回のハミルトニアン行列の対角化で求めることができる。また、SAC/SAC-CI 法は
size-extensive/size-intensive であり[4]、大きな共役系の計算にも適用可能である。
Fig. 1 に SAC-CI による計算結果と実験の UV-vis スペクトルを示す。対称性から、ほとんどの状態が光学禁制であるが、実験では、振電相互作用による弱い吸収が可視光領域に観測されている。可視光領域には、11T1g、11T2g、11Gg、11Hg 状態が計算され、これは、縮退した HOMO→LUMO
遷移に相当する。Fig. 2 のように C60 の HOMO 6hu は 5 重に、LUMO7t1u は 3 重に縮退しており、HOMO→LUMO 遷移は、T1g + T2g + Gg + Hg の 15 状態に相当する。軌道レベルで縮退した状態は、配置換相互作用により、そのエネルギー準位が変わる。1 電子励起 CI の準位は 11T2g < 11T1g < 11Gg
< 11Hg であった。SAC-CI では、さらに電子相関が考慮される。縮重度の大きな状態ほど電子相
関による安定化が大きく,SAC-CI による準位は 11Gg < 11T2g < 11T1g < 11Hg であった。
400 nm (3.10 eV)付近に band A と呼ばれる弱い吸収が観測される。従来の計算では、これは、光学許容な 11T1u に帰属されてきた。しかし、我々の計算では、これは、光学禁制の 11T2u に帰属される。これは、HOMO→LUMO+1 遷移で、弱い吸収強度は、振電相互作用によるものと推測される。SAC-CI の結果では、11T2u は、他の HOMO→LUMO+1 遷移である 11Gu や 11Hu とエネルギ
ー的に十分離れており、孤立したピーク(band A)として観測されることと矛盾しない。また、光学
許容の 11T2u は band B に、21T2u は band C に帰属することで、実験スペクトルを矛盾無く説明できた。また、11Gu や 11Hu 等のその他の禁制状態の振動吸収は、11T2u とともに band B として観測されていると考えられる。
SAC-CI 計算に基づいた、スペクトルの帰属は、従来の半経験的計算及び pure GGA 汎関数による TDDFT 計算に基づいた帰属とは一部異なる。SAC-CI の結果は、むしろ PBE0 のような hybrid
汎関数による TDDFT の結果に近かった。また、複雑に縮退した電子状態の信頼できる記述には、
電子相関を適切に考慮することが必要である。

Fig. 1. Experimental absorption spectrum of C60 in
n-heptane [5] and the SAC-CI computational results.
The vertical lines represent optically allowed states.
The cross symbols and open circles denote opticallyforbidden gerade (even) and ungerage (odd) states,
respectively.

Fig. 2. Valence orbital structure of C60 in Ihsymmetry.

[1] G. Orlandi and F. Negri, Photochem. Photobiol. Sci. 1, 289 (2002).
[2] R. Fukuda and H. Nakatsuji, J. Chem. Phys. 128 094105 (2008).
[3] R. Fukuda and M. Ehara, J. Chem. Phys. 136, 114304 (2012); R. Fukuda, M. Ehara, and H. Nakatsuji, J.
Chem. Phys. 133, 144316 (2010).
[4] H. Nakatsuji, T. Miyahara, R. Fukuda, J. Chem. Phys. 126, 084104 (2007).
[5] J. Catalán, Chem. Phys. Lett. 223, 159 (1994).

3E12

凝集系における分子の励起状態：
局在化軌道を用いた分子間相互作用の解析(京大福井セ 1, 神戸大情報 2)

長谷川 淳也 1, 藤本

和宏 2, 河津励1

A LMO picture for excited-state molecular interactions inproteins and solutions(Kyoto Univ.1. Kobe Univ.2) J. Hasegawa1, K. J. Fujimoto2, and T. Kawatsu1【序】周辺環境に応じて分子の励起エネルギー準位が変化する現象はソルバトクロミズムとして良く研究されてきた。我々はレチナール蛋白質、蛍光蛋白質、蛍ルシフェリンなどの光機能性蛋
白質におけるカラー・チューニングの起源について、QM(SAC-CI)/MM(AMBER)計算に基づく研究を行い、周辺環境の静電ポテンシャルが重要な役割を果たすことを明らかにした[1]が、他方で、周辺環境の電子状態が無視できない影響を与えることも見出している[1,2]。バクテリオロドプシン(bR)及びメタノール溶液中におけるレチナール色素の第一励起状態について、約 6 Å 程度の距
離にあるアミノ酸までを考慮した CIS 計算を行うと、周辺環境を点電荷モデルで記述した結果と比較して、bR では約 0.3 eV、MeOH では約 0.1 eV 励起エネルギーが減少する。この差は、ソルバ
トクロミズムの物理化学的原理を議論する上では無視できない寄与を与える[2]。
我々は周辺環境の電子状態に由来する効果を解析するために局在化分子軌道を計算する手法を提案した[3]。一電子軌道をフラグメント（アミノ酸残基やペプチド結合）へと帰属できるので、励起状態の波動関数の解析が明確になった。本研究では、同局在化軌道を用いた解析を行い、凝集系における分子の励起状態の波動関数の構造を明らかにすることを目指す。
【計算】周辺環境の電子状態がレチナール色素の励起エネルギーに及ぼす効果を記述するために、溶媒分子（以降、蛋白質も含めて周辺環境を溶媒と総称する。）を含めた分子クラスター（図 1）を構築し、CIS 計算を行った。クラスターの外側の溶媒分子の効果は静電ポテンシャルとして考慮した。このモデルでは、溶媒―色素間の電荷移動（CT）励起（軌道間相互作用）、溶媒分子における局所励起（分子内分極）、溶媒分子間の CT（分子間分極）が考慮される。分子クラスターのサイズ拡大に伴い第一励起エネルギーは減少するが、およそ 6 Å 半径（約 1000 原子系）で、励起エネルギーの減少は収束した。分子軌道の局在化は、我々が提案した方法[3]を用い、この一電子基底を用いて CIS 計算を行なった。実験結果と比較する際は、ONIOM 法を用いて、SAC-CI 法
により計算した色素の励起エネルギーに CIS 法による環境効果を外挿した結果を用いた。
【結果と考察】周辺環境の効果が最も顕著であった bR について、局在化軌道を用いた解析結果について説明する。図 1 には第一励起状態の主配置に関与する分子軌道の一部を示す。最も主要な
配置はレチナール色素の HOMO（図 1a）から LUMO（図 1b）への一電子遷移であり、これに Tyr185の HOMO（図 1c）からレチナールの LUMO への CT 遷移が副配置として混合することが明らかになった。この結果は、レチナールの HOMO が Tyr の HOMO と相互作用し、非局在化すること

により励起エネルギーが低下したことを意味する。
次に、より広域な電子状態の変化を理解するために、電子励起をフラグメント単位に分類し、(a) MO1746: HOMO of PSB(b) MO5014: LUMO of PSB(c) MO1491: HOMO of Tyr185

図 1. bR の第一励起状態の波動関数に関与する局在化分子軌道。(a)Protonated retinal Schiff base (PSB)の
HOMO 及び (b) LUMO。(c) Tyr185 の HOMO。

各グループにおいて最大の CI 係数を抽出した。結果を図 2 に示す。bR では主配置である PSB における励起に加えて、Tyr185 から PSB への CT（及び逆 CT）励起、Tyr185 における局所励起が比較的大きな係数を持っている。他方で、MeOH 溶液中の場合は、レチナールから溶媒分子への CT
励起が見られるが、bR ほど顕著な寄与が無いことが特徴的である。
(a)

(b)

図 2. (a)bR 及び(b)MeOH 中のレチナール色素について、第一励起状態の波動関数をアミノ酸残基毎に分類し、最大の係数をプロットした。

また、図 2a と 2b では、電子励起の寄与がレチナールから溶媒分子への CT 励起、その逆の CT
励起、エキシトン励起を表わす３つの直線上に偏って分布することが明確に分かる。この状況は他の系でも共通であり、周辺環境の励起波動関数の構造を示唆している。
図 3 には、各フラグメントについて励起波動関数のノルムを計算し、PSB からの距離に関して、プロットしたものを示す。bR において、距離約 4 Å 程度の第一溶媒和圏では CT 等の効果が顕著に見られ、軌道間相互作用の重要さが示唆された。それを超える領域では CT の寄与は減衰し、エキシトンの寄与が相対的に重要になる。遷移
双極子モーメント間相互作用の長距離性が反映していると考えられる。

【参考文献】[1] J. Hasegawa, K. J. Fujimoto, and H. Nakatsuji
ChemPhysChem 12, 3106 (2011).
[2] K. Fujimoto, K. Asai, and J. Hasegawa, Phys.
Chem. Chem. Phys. 12, 13107 (2010).
[3] J. Hasegawa, T. Kawatsu, K. Toyoya, and K.
Matsuda, Chem. Phys. Letters, 508, 171 (2011).

図 3. 各フラグメントの局所励起配置、色素とのCT 励起配置のノルムを、フラグメント―色素間距離についてプロットした。（右）第一溶媒和圏
において主要な寄与を与えるフラグメント。

3E13
高信頼性計算分子科学技術－キラサク(吉良作)－による
円二色性スペクトルの構造依存性に関する研究（量子化学研究協会 1, JST-CREST2）宮原友夫 1,2, 中辻博 1,2

Conformational dependence of circular dichroism spectra: CHIRASAC study（QCRI1, JST-CREST2）Tomoo Miyahara1,2, Hiroshi Nakatsuji1,2【序】円二色性（CD）スペクトルは、溶液中の構造を調べるのによく利用されている。しかし、その解釈は定性的であり、問題も多い。従って、キラル分子や、キラル分子と相互作用したアキラル分子の円二色性
スペクトルの理論計算・解析は、構造を解明・予測するために重要である。これまでの研究から、SAC-CI
法[1-3]が円二色性スペクトルの研究に有用であり、また、そのスペクトルは構造変化に極めて敏感である
ことが分かっている。そこで、キラル分子の構造解析と機能創出を目的として、DFT 法の構造計算とSAC-CI 法による円二色性スペクトル計算を中心とした、高信頼性計算分子科学技術―キラサク（吉良作）－を構築することを目指している。（図１）【方法】Gaussian[4]に搭載されている DFT 法と SAC-CI 法を用いた。下で述べるキラル分子を DFT法により構造最適化し、その構造を用いて SAC-CI 法による円二色性スペクトル計算し、実験スペクトルと比較検討することにより、キラル分子の溶液中の安定構造、円二色性スペクトルの立体配座依存性、相互作用情報などについて考察する。
【Dihydroazulene 誘導体】は 6 位にメチル基がある誘導体(A)と 6 位が水素の誘導体(B)の実験吸
収(UV)・円二色性(CD)スペクトルが報告されている[5]。
実験吸収(UV)スペクトル(図 2(a,c)黒)は、誘導体 A と B でほとんど変化がないが、実験円二色性(CD)スペクトル(図 2(b,d)黒)は誘導体 A では 240nm に小さな正のピークが現れるのに対して、誘導体 B ではピークが現れない。従って、

図1 高信頼性計算分子科学技術「キラサク（吉良作）」

図 2. Dihydroazulene 誘導体(A,B)の SAC-CI UV・
CD スペクトル(赤)と実験 UV・CD スペクトル(黒)[5]。

円二色性スペクトルから置換基の有無を判別することができる。他方 SAC-CI UV・CD スペクトルは、共に実験スペクトルを良く再現している(図 2 赤)。SAC-CI UV スペクトルでは、実験で観測されたエネルギー範囲より高い励起状態でスペクトルに違いがあることを示唆している。一方、
SAC-CI CD スペクトルでは、240 nm のピークに対応する 2 つの励起状態が計算された。誘導体
A では共に正の値を持つのに対し、誘導体 B では一方が負の値、もう一方は強度が弱く、結果として、誘導体 A のみ正のピークが観測されるということが明らかになった。また、320nm 付近の
正のピークと 220 nm 付近の負のピークで、SAC-CI CD スペクトルは、誘導体 A と B の強度比も正しく再現している。このように SAC-CI 法により円二色性スペクトルを本質的に理解することによって、分子構造や相互作用の nature を知ることができるのである。
【N-(1-phenylethyl)thiazoline-2-thione】の実験円二色性スペクトル[6]は温度に依存し、高温(+22 度)と低温(−148 度)で強度が異なり、低温では 350nm 付近に負の符号のピークが現れる(図3(b,c))。気相中の最安定構造を元にして、図3(a)の矢印で示すように、キラル炭素と結合したフェニル基とthiazoline をそれぞれ
30 度ずつ回転させた
SAC-CI CD スペクトルを計算した。最安定構造の SAC-CI CD ス
ペクトルは 2 個の実験
ス ペク トル のど ちらと も一 致し なか った(図 3(d))。一方、構造A と B の SAC-CI CD
スペクトル(図 3(e,f))
は、それぞれ高温と低図 3. (a) N-(1-phenylethyl)thiazoline-2-thione、(b)-(c)実験円二色性(CD)スペクトル[6]、(b) +22 度、(c)−156 度、(d)-(f) SAC-CI CD スペクトル、(d) 最安定構造：φ1 = −165.0, φ2 = +46.5、(e) φ1 = −165.0, φ2 = +166.5、(f) φ1 = +165.0, φ2 = +196.5。

温の実験スペクトルと一致した。このことから、高温では構造 A、低温では構造 B のような立体配座をしていることが推測できる。実験スペクトルは混合溶液(diethyl ether, isopentane,ehanol)中で測定されているが、このような混合溶液の溶媒効果を計算に簡易的に含めることは困難である。しかし、SAC-CI CD スペクトルの計算から、溶液中の立体配座を予測することができる。詳細は会場で述べる。
【参考文献】[1] H. Nakatsuji, Chem. Phys. Lett.1978, 59, 362.; 1979, 67, 329, 334; Bull. Chem. Soc. Jap. 2005, 78, 1705.
[2] M. Ehara, J. Hasegawa, H. Nakatsuji, Theory and applications of Computational Chemistry, The First 40
Years, Elsevier Oxford, 2005; p1099.
[3] SAC-CI homepage. http://www.qcri.or.jp/sacci/ (6/6/2005)[4] M. J. Frisch, et at. GAUSSIAN 03, Gaussian, Inc. Pittsburgh PA, 2003.
[5] N. Harada, J. Kohori, H. Uda, K. Nakanishi, R. Takeda, J. Am. Chem. Soc. 1985, 107, 423.
[6] J. Roschester, U. Berg, M. Pierrot, J. Sandström, J. Am. Chem. Soc. 1987, 109, 492.

3E14

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

An Improvement in Solute-Solvent Interactions via the
Number-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 が以下のように与えられる。
M

NT 1
d 2r
 a  r  t   V a  r  t  



dt 2
a 1

 NT 1 a    r  t    1
 a1


NT 1

NT 1

V ad  r  t      a  r  t  V a  r  t       a  r  t  V a  r  t  dr
a1

r

a1

(1)
(2)

4

4

SD (OC-Ow )

3.5
3

2

1.5

1

1

0.5

0.5

Cα

Cc

2

1.5

0

OC2

NAM (N-Hw )

2.5

g(r)

g(r)

SD (N-Hw )NAM (N-Ow )3

NAM (OC-Hw )2.5

SD (N-Ow )3.5

SD (OC-Hw )NAM (OC-Ow )

Hα
N

O1C
H 3N

H1,2
N

0
1

2

3

4

5

6

7

r (Å)
(a) -COO 基(OC)-水分子(OW,HW)間

1

2

3

4

5

6

7

r (Å)

(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).