1D15

MF-QM/MM法による化学反応の自由エネルギー計算と溶液・酵素反応への応用（京大院・理）山本 武志，中農 浩史
Free energy calculation of chemical reactions using MF-QM/MM method
and its application to solution and enzyme reactions(Kyoto University) Takeshi Yamamoto, Hiroshi Nakano【序】計算機性能の進歩と電子状態理論の発展によって、様々な分子のエネルギーを正確に求める事が可能になってきている。特に、限られた数の構造サンプルに対してエネルf the QM molecule
given by
ギー計算を行う場合、非常に高精度な電子状態理論を用いる事が出来る。一方、系の時間発展を追跡したい場合や、自由エネルギーなど統計平均を必要とする場合、電子状態計算
nucele!
!
Zδ(x−R)−
δ(x − xi ),(3)ρ̂(x)=の回数が非常に多くなるため、コストの高い理論を用いるのは事実上困難になる。すなわ
aacule given byai
ち、ダイナミクスや自由エネルギーの計算では、計算レベルと計算コストに関するジレン
nucele!
s the electrostatic
potential!
(ESP) generated by the MM molecules, namely,
マの問題が大きく、それを軽減するための方法が模索されている。同じ問題はいわゆるZa δ(x − Ra ) −δ(x − xi ),(3)
ρ̂(x) =a QM/MM計算（系の一部を電子状態理論で記述し、残りを経験的力場で記述する方法）でiMM! qi.
(4)v(x)=も同様に起きる。特に、化学反応のQM/MM自由エネルギー計算は負荷が大きいため、こatic potential(ESP) generated by the|xMM
molecules,
namely,−r|ii
基盤Ｃ（一般）‐２
れを現実的なコストで行うための方法が研究されている。
MM!
qi 目 involves
calculation of the aboveA(R)repeated
研 究的（つづき ）  QM calculation of Ψ for a large.
(4)v(x) =|x − ri | (e.g., more than 104 ). To reduce the computationalQM/MM 法で計算したいとションでも発生する。例えば、溶液反応の自由エネルギーを
different
solvent
configurationsi【方法】このような背景から、発表者らはQM波動関数に統計的な平均場近似を入れるタ
perator of the QM molecule given by
すると、溶媒の様々な配置に対して溶質の電子状態計算を繰り返し行わなければならず、イプのアプローチ(ここではMF-QM/MMと呼ぶ）を研究している。QM/MM計算に部分
onsider
replacing
the計算コストが非常に大きい。この問題を緩和するために、申請者はこれまで自由エネルギ
exactwave!
functionΨ(R, r) by
one, Ψ̃(R), which
he aboveA(R)involves
repeated
QM calculationofΨ fortrial
a largenucelesome!
的な平均場近似を入れるアイデアは昔からいくつかのグループによっ− theδ(x− xone.
(3)
ρ̂(x)
= 104Za δ(x − Ra ) toi ),ーの変分原理を用いて、QM/MM
自由エネルギーを効率よaverageor(e.g.,coarse-grained
approximationexact
One can thenntsomeconfigurations
morethan). To reduce thecomputational
自由エネルギープロファイル
a の 方 法 論 を 開 発iし て き た （ Yamamoto,
く求めるためて使われているが([1,2]のIntroductionを参照）、その妥当性
approximate
total functionenergy inΨ(R,termsasng
the exact waver) ofbyΨ̃(R)sometrial one, Ψ̃(R), whichJ.Chem.Phys.2008）。この方法では、溶質の電子状態計算
nd v(x) is theelectrostatic
potential(ESP)
generated by the MM molecules, namely, 厳密なQM/MM 計算
や根拠はあまりはっきりしておらず、精度に関する検証も不十"を溶媒の各配置について行う代わりに、溶媒が作る平均的近似的な QM/MM 計算（本研究）
or coarse-grained approximation
to the exact one. vdW
One can thenMMẼtot (R,分であった。そこで発表者らはまず自由エネルギーの変分法r; Ψ̃) = な場のもとで溶質の電子状態計算を行う。また、これに対"Ψ̃|Ĥ0 + dxρ̂(x)v(x)|Ψ̃# + EQM/MM + EMM ,(5)!
RISM-SCF 法
qital energy interms of Ψ̃(R) as.
(4)v(x) =を用いてQM/MM自由エネルギーと平均場近似の関係を定式応して溶媒の分布を再計算し、溶質の電子状態と溶媒の分誘電体モデル(PCM,COSMO)
|x − ri |"irresponding化した
approximate
free energy as[1]。具体的には次のような量を考える：
布がセルフコンシステントになるようにする。この方法でvdWΨ̃) = "Ψ̃|Ĥ0 + dxρ̂(x)v(x)|Ψ̃#
EQM/MM+E ,(5)" +A(R)は、電子状態計算の回数はたかだか数十回で済むため、厳Rigorous calculationof the1above
involves MM
repeated QM calculation of Ψ for a largeNln drサンプリング
exp[−β Ẽtot(R, r; Ψ̃)].
(6)
Ã(R; Ψ̃)密な= −QM/MM（数万回の電子状態計算を必要
βumberofdifferent
solvent
configurations(e.g., more than 104 ). To reduce the computational
pproximate free energyasとする）に比べると、計算コストが格段に小さくなる。ま
"
ost,
we
considerreplacingwave function
Ψ(R,by sometrial one,Ψ̃(R), whichÃ above
dependparametricallyon Ψ̃.
To r)proceedfurther,
we noteboth Ẽtot and the exactた、ベンチマーク系を用いたテスト計算から、得られる自1ここでA(R;N Ψ) は溶質が構造Rを持つときの全系のHelmholtzexp[−β Ẽtot (R, r; Ψ̃)].
(6)R; Ψ̃) = − ln dr由エネルギー曲面は厳密な
QM/MM 計算の結果と非常に
epresentssome
β average or coarse-grained approximation to the exact one. One can thenlity
よく一致することを示すことが出来る (Nakano et al, 投
自由エネルギーで、E(tot)は系の全エネルギーである(rはMM
efineanapproximate
total
energyterms
Ψ̃(R)asr)on≤inẼtot(R,proceedr;ofΨ̃),
Etot (R,稿中)。概念的には、上の方法は厳密なQM/MMÃ above depend
parametricallyΨ̃.
Tofurther,
wenote 法と従来 (7)部分の座標）。E(tot)として次のような形を仮定する："の近似的な方法（誘電体モデルや RISM-SCF 法など）の中間的な方法と言える（図・上）。
vdWds true for arbitraryΨ̃(atanygivenRandr)becauseΨ̃isnot
necessarilytheẼtot (R,r; Ψ̃) = "Ψ̃|Ĥ0 + dxρ̂(x)v(x)|Ψ̃# + EQM/MM + EMM ,(5)必要な電子状態計算の回数が少ないことから、より信頼性の高い電子状態理論（MP2法#やCCSD(T)法など）を用いて自由エネルギー計算を行うことが出来る（図・下、典型的(R, r) ≤Hamiltonian,Ẽtot (R, r; Ψ̃),Ĥ0 + dxρ̂(x)v(x). The above(7) relation thenEtoteffective
ate of anな溶液内SN2反応）。このような特性を生かして、現在では、従来計算が困難であったnd the corresponding
approximate
freeenergyas
ここで、Ψ は溶質の（統計的な）試行波動関数、v(x)はMM部分がQM部分にかける（瞬trary Ψ̃ (at any givenand r) because
ATPR分子を含む溶液中リン酸エステルの反応自由エネルギー計算を行なっている。
" Ψ̃ is not necessarily the# A(R) ≤ 1Ã(R; Ψ̃), Nを系の厳密なQM/MM自由エネルギーとする間的な）静電ポテンシャルである。A(R)(8)
dr exp[−βẼtot (R, r;thenΨ̃)].
(6)= − ln The
tive Hamiltonian, Ĥ0Ã(R;+ Ψ̃)dxρ̂(x)v(x).
above relationβ
本研究の目的：励起状態ダイナミクスのための新規なQM/MM計算手法の開発と応用
 は次の不等式を満たすことが分かる[1,2]。
と、 A(R;Ψ)
es that Ã(R; Ψ̃) gives an上の【背景upper boundon the exact A(R) irrespective of the choice
A】で述べたように、凝縮相中の励起状態ダイナミクスをQM/MM 法で取Note that both Ẽtot and Ã above depend parametrically on Ψ̃. To proceed further, we note
り扱うには２つの問題点（電子状態理論の信頼性と統計サンプリングの問題）が存在する。
A(R)Ã(R; Ψ̃), best trial wavefunction Ψ̃ by minimizing(8)
can then obtainthe≤statisticallyÃ(R; Ψ̃)he inequality
そこで本研究では、必要な電子状態計算の回数を大幅に削減するための新規な QM/MMcttoΨ̃.
Theresulting
freeenergy
isA(R)
writtenas
計算手法を開発し、それを通じて信頼性の高い電子状態理論の使用を可能にすることを目) gives an upper boundon theexactirrespectiver) ≤Ẽtot (R, r; ofΨ̃),the choice(7)
Etot (R,的とする。
（これは【背景 B】で述べたこれまでの研究の一種の拡張に相等する）。そうし

n the statistically best trial
Ψ̃ byminimizing Ã(R; Ψ̃)
AMFwavefunction
Ã(R;Ψ̃).
(9)
(R)
min
参照結果と比較することで方法の妥当性を検
which holds true for て得られた方法をベンチマーク系に適用し、arbitraryΨ̃ =(atanygiven
R and r) becauseΨ̃ is not necessarilytheΨ̃#esulting free energy 証する。また、溶液内包接色素分子へ応用し、空間的に制限された色素分子の励起状態ダis written as
round state of an effective Hamiltonian, Ĥ0 + dxρ̂(x)v(x).
The above relation then˜n then obtain the statistically best trial wavefunction Ψ̃ by minimizing Ã(R; Ψ̃)ling function for the MM part defined asto Ψ̃. The resulting free
" energy is written as

vdW+ EMM ,(12)EMM (R, r; ρ̃) = dxρ̃(x)v(x) + EQM/MMAMF (R) = min Ã(R; Ψ̃).
(9)
 をΨ Ψ̃について最小化すればベストな自由エネルギーの近似が得られるこΨ)
従って、 A(R;he charge density derived from Ψ̃, namely, はMM部分が作る平均的な場の下で計算することになる。
ng the variational
procedure, we obtain
the equation for Ψ̃ as follows:とが分かる。また、Ψ$%"ρ̃(x) = !Ψ̃|ρ̂(x)|Ψ̃".
(13)
Ĥ0 + dxρ̂(x)"v(x)# |Ψ̃# = EQM |Ψ̃#.
(10)ng Ψ̃ thus obtained into AMF (R) in Eq. (9), we obtain
その結果、ベストな近似的自由エネルギーは下の形で得られる[1,2,4]：5(14)AMF (R) = !Ψ̃|Ĥ0 |Ψ̃" + ∆µ(R, ρ̃),

上の方法は、明らかに従来の誘電体モデルやRISM-SCF法をQM/MM法に拡張したものに
sical contribution to the solvation free energy ∆µ is given by?
なっている。この方法の利点は以下の通りである：（１）誘電体モデルと同程度のQM計"
1算で自由エネルギーが求められるため、比較的高精度な電子状態理論を用いることが出来∆µ(R, ρ̃) = − ln drN exp[−βEMM (R, r; ρ̃)].
(15)
βる；（２）QM計算のステップとMMサンプリングのステップをself-consistentになるまで繰り返す必要があるが、MMサンプリングのステップは完全に古典的な計算となるた
ional simplicity,
we neglect the partition function of the bulk MM system
Δµ (R, ρ ) の評価
tion above.)め、高速なMDプログラムやハードウェアを利用する事が出来る；（３）Obviously, the above expression for A (R) is analogous toMF

をMMモデルで行うため、均一な溶液系だけでなく、酵素など特異な反応部位を持つ系に
raditional solvation
models (e.g., the PCM and RISM-SCF methods). The
ence is that も使うことが出来る[5,6]。また、溶媒だけでなくタンパク質の揺らぎを考慮出来る。
one samples solvent configurations explicitly, e.g., by performingn in the presence of solute charge density ρ̃(x). Note that the solute charge【結果】ここでは、まず上の方法をいくつかの典型的な溶液反応に適用し、厳密なQM/
fixed during the MD calculation, so that that the ensemble average in Eq. (11)MM計算と比較することでその信頼性について検証した。単純なSN2反応では、MF-QM/
assical in thesense that no repeated QM calculations are involved. We note
MM計算によって求めた自由エネルギープロファイルは厳密なQM/MM結果と非常に良く まわりの）統計揺らぎが自由エネルeffective Schrödinger
equation in Eq. (10) is nonlinear with respectto Ψ̃,一致することが分かる。これは、真の波動関数の（Ψdepends onギーにほとんど影響しないことを示している（詳細はポスター発表[7]を参照）。次に、
ρ̃(x) = !Ψ̃|ρ̂(x)|Ψ̃" via the ensemble average in Eq. (11). This

より現実的な反応として、水中のphosphate
diesterfor
monoanionとその誘導体の加水分
e needs to performthe QM calculation of Ψ̃ and the MD sampling!v(x)"解 [4]、および、PchB中のisochorismate
pyruvate lyase酵素反応への応用を行い、そe manner until
the self-consistent Ψ̃ is obtained.の結果について説明する（左図：phenyl methyl phosphateのMP2レベルでのPMF、中・右図：isochorismateとPchB）

sion to periodic QM/MM systems(b) 30(a) 3020
we have assumed that the system consists of one
solute molecule and N
PMF (kcal/mol)

PMF (kcal/mol)

201010les under nonperiodic
boundary condition. Here, weconsider how to modify
00−10−20−30

6

MF−QM/MM(dens)
MF−QM/MM(sc)MF−QM/MM(pc) N.A.
−2−101! = r(P−Olg) − r(P−Onuc) (Å)−10−20−30

2

MF−QM/MM(dens)
MF−QM/MM(sc)MF−QM/MM(pc) N.A.
−2−101! = r(P−Olg) − r(P−Onuc) (Å)

2

9: PMF profilesJ.Chem.Phys.
of reaction of hydroxide
ion244104
attack on methyl
phenyl phosphate monoanions[1] FIG.
T. Yamamoto,129,(2008)
solution.
They are obtained
from the Phys.
mean-field
QM/MMcalculations
with charge[2] inH.aqueous
Nakano,T. Yamamoto,J. Chem.
136,134107(2012)(dens) andscreened charges(sc). ThePMFLett.
cannot(2012,
be obtained
the calculation[3] density
H. Nakano,T.withYamamoto,
Chem.
Phys.
in frompress)[4] with
H. Nakano,T.because
Yamamoto,
J. Chem.
TheoryComput.
(submitted)
point chargesof their divergent
behaviour.
The QM
calculationis performed at the[5] MP2/6-311+G(2d,2p)
T. Kosugi, S. Hayashi,
J. Chem. Theory Comput. (2011)
level.
[6] T. Kosugi, S. Hayashi, J. Am. Chem. Soc. 134, 7045 (2012)[7] 1P097 QM/MM計算における波動関数の統計揺らぎの効果と平均場近似の妥当性の検証