QM/MM Free Energy Perturbation (FEP) 法：基礎編

化学反応の解析において、量子力学 (QM) と分子力学 (MM) を組み合わせた QM/MM 法は、計算の効率性と精度を両立させるパワフルな手法として活用されています。

反応経路に沿った自由エネルギー変化を計算する際に、私は QM/MM 法と Free Energy Perturbation (FEP) 法を組み合わせた手法をよく用います（Yamamoto 2020, 2021, 2022）。QM/MM FEP 法は、系のエネルギー状態の変化を小さな摂動として扱い、その変化を統計力学的に解析することで、反応経路に沿った自由エネルギー変化を求める手法です。

このノートでは、私がよく使っている QM/MM FEP 法について、ざっくりと解説します。

QM/MM 法の基本構造

QM/MM 法の起源は1976年の Warshel と Levitt の論文にさかのぼります（Warshel 1976）。彼らは酵素反応の計算に初めてこの手法を適用し、Nobel 賞につながる成果を挙げました。以来、計算機の発展とともに手法も進化し、今日では様々な分野で不可欠なツールとなっています。

nobelprize.org
nobelprize.org

QM/MM のエネルギー分割

QM/MM 法では、系全体のエネルギー EtotalE_{\text{total}} は、QM 部分、MM 部分、そしてその間の相互作用エネルギーから成り立っています：

Etotal=EQM(rQM)+EMM(rMM)+EQM/MM(rQM,rMM) E _\text{total} = E _{\text{QM}}(\mathbf{r} _\text{QM}) + E _\text{MM}(\mathbf{r} _\text{MM}) + E _\text{QM/MM}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{MM})

  • EQM(rQM)E_{\text{QM}}(\mathbf{r}_\text{QM}) ：QM 部分（反応に直接関与する部分）のエネルギーであり、通常、密度汎関数理論（DFT）法や Hartree-Fock 法などを用いて計算します
  • EMM(rMM)E_{\text{MM}}(\mathbf{r}_\text{MM}) ：MM 部分（環境や溶媒など）のエネルギーであり、分子力学ポテンシャルに基づいて算出します
  • EQM/MM(rQM,rQM)E_\text{QM/MM}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM}) ：QM 部分と MM 部分間の相互作用エネルギーであり、電荷相互作用やファンデルワールス相互作用が含まれます

QM/MM 間の相互作用

QM 部分と MM 部分の相互作用は、以下のように分割できます：

EQM/MM=Ees(rQM,rQM)+EvdW(rQM,rQM)+EMM-bonded(rQM,rQM) E _{\text{QM/MM}} = E _{\text{es}}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM}) + E _{\text{vdW}}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM}) + E _{\text{MM-bonded}}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM})

  • Ees(rQM,rQM)E_{\text{es}}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM}) ： QM 部分と MM 部分の間の静電相互作用
  • EvdW(rQM,rQM)E_{\text{vdW}}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM}) ： ファンデルワールス相互作用
  • EMM-bonded(rQM,rQM)E_{\text{MM-bonded}}(\mathbf{r} _\text{QM}, \mathbf{r} _\text{QM}) ： MM 部分と QM 部分にまたがる結合に関連するエネルギー

Free Energy Perturbation 法

QM/MM Free Energy Perturbation (FEP) 法は、Zhang らによって効率的な反応経路の自由エネルギー変化を評価するために導入されました（Zhang 2000）。この手法は、反応座標に沿った QM 部分と MM 部分の相互作用を考慮しながら、統計力学的に正確な自由エネルギー差を得ることを可能にします。彼らのアプローチでは、QM 部分を反応座標に固定し、MM 部分の揺らぎを統合することで、QM/MM 系全体の自由エネルギーを求めています。

配置分配関数と自由エネルギー

化学反応において、反応経路に沿った自由エネルギー変化を理解するためには、系の分配関数を計算することが重要です。系の配置分配関数 ZZ は以下のように定義することができます：

Z=∫exp⁡(−β[EQM(rQM)+EMM(rMM)+EQM/MM(rQM,rMM)])drQMdrMM Z = \int \exp\left( -\beta \left[ E _{\text{QM}}(\mathbf{r} _{\text{QM}}) + E _{\text{MM}}(\mathbf{r} _{\text{MM}}) + E _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}, \mathbf{r} _{\text{MM}}) \right] \right) d\mathbf{r} _{\text{QM}} d\mathbf{r} _{\text{MM}}

ここで、β=1kBT\beta = \frac{1}{k_B T} は逆温度、rQM\mathbf{r} _{\text{QM}} と rMM\mathbf{r} _{\text{MM}} は、それぞれ、QM 部分と MM 部分の座標を表しています。

反応座標 RcR_c に関連する部分分配関数 Z(Rc)Z(R_c) は次のように表すことができます：

Z(Rc)=∫δ(Rc−f(rQM))exp⁡(−β[EQM(rQM)+EMM(rMM)+EQM/MM(rQM,rMM)])drQMdrMM Z(R _c) = \int \delta(R _c - f(\mathbf{r} _{\text{QM}})) \exp\left( -\beta \left[ E _{\text{QM}}(\mathbf{r} _{\text{QM}}) + E _{\text{MM}}(\mathbf{r} _{\text{MM}}) + E _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}, \mathbf{r} _{\text{MM}}) \right] \right) d\mathbf{r} _{\text{QM}} d\mathbf{r} _{\text{MM}}

この分配関数に基づく自由エネルギー F(Rc)F(R_c) は次のように定義できます：

F(Rc)=−1βln⁡Z(Rc) F(R_c) = -\frac{1}{\beta} \ln Z(R_c)

QM と MM の分離

QM 部分と MM 部分の自由度を分離して扱うために、MM 部分に関する自由エネルギーの平均力ポテンシャルを導入します。これを式で表すと以下のようになります：

FQM/MM(rQM)=−1βln⁡∫exp⁡(−β[EMM(rMM)+EQM/MM(rQM,rMM)])drMM F _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}) = -\frac{1}{\beta} \ln \int \exp\left( -\beta \left[ E _{\text{MM}}(\mathbf{r} _{\text{MM}}) + E _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}, \mathbf{r} _{\text{MM}}) \right] \right) d\mathbf{r} _{\text{MM}}

この式では、MM 部分の統計力学的な揺らぎを取り込んでいます。したがって、反応座標に沿った自由エネルギープロファイルは以下のようになります：

F(Rc)=EQM(rQMmin)+FQM/MM(rQMmin)−1βln⁡∫δ(Rc−f(rQM))exp⁡(−βϵ(rQM))drQM F(R _c) = E _{\text{QM}}(\mathbf{r} _{\text{QM}}^{\text{min}}) + F _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}^{\text{min}}) - \frac{1}{\beta} \ln \int \delta(R _c - f(\mathbf{r} _{\text{QM}})) \exp\left( -\beta \epsilon(\mathbf{r} _{\text{QM}}) \right) d\mathbf{r} _{\text{QM}}

ここで、rQMmin\mathbf{r} _{\text{QM}}^{\text{min}} は反応座標に沿った最適な QM 部分の座標であり、ϵ(rQM)\epsilon(\mathbf{r} _{\text{QM}}) は QM 部分の振動運動に対応するエネルギーの変動です。

FEP 法による自由エネルギー変化

反応経路に沿った自由エネルギーの相対変化は、QM 部分の振動が反応座標に沿って一様に変化すると仮定すると、以下のように求められます：

ΔF(Rc)≈ΔEQM(rQMmin)+ΔFQM/MM(RcA→RcB) \Delta F(R_c) \approx \Delta E _{\text{QM}}(\mathbf{r} _{\text{QM}}^{\text{min}}) + \Delta F _{\text{QM/MM}}(R_c^A \to R_c^B)

ここで、ΔFQM/MM(RcA→RcB)\Delta F_{\text{QM/MM}}(R_c^A \to R_c^B) は次のように定義されます：

ΔFQM/MM(RcA→RcB)=−1βln⁡⟨exp⁡[−β(EQM/MM(rQMmin(RcB))−EQM/MM(rQMmin(RcA)))]⟩MM,A \Delta F _{\text{QM/MM}}(R_c^A \to R_c^B) = -\frac{1}{\beta} \ln \left\langle \exp \left[ -\beta (E _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}^{\text{min}}(R_c^B)) - E _{\text{QM/MM}}(\mathbf{r} _{\text{QM}}^{\text{min}}(R_c^A))) \right] \right\rangle _{\text{MM}, A}

ここで、RcAR_c^A と RcBR_c^B は反応座標の異なる点を表し、⟨⋯ ⟩MM,A\langle \cdots \rangle_{\text{MM}, A} は QM 部分を RcAR_c^A に固定した状態での MM 部分に関する平均です。この式から、QM 部分と MM 部分の相互作用エネルギーの変動に基づいて、自由エネルギーの変化を計算することができます。

References

Posted :