メインコンテンツへスキップ

メタダイナミクスのパラメータ設定のコツ

  • 更新

メタダイナミクス法（Metadynamics）は、自由エネルギー地形の極小点に留まってしまう系に対し、ガウス型のバイアスポテンシャルを積み上げることで強制的にエネルギー障壁を乗り越えさせる強力な手法です。

しかし、効率的かつ正確な計算を行うためには、ガウス関数の高さや幅、バイアスファクターなどのパラメータを適切に設定する必要があります。

本記事では、PLUMEDを用いたWell-Tempered Metadynamicsを実行する際のパラメータ設定の経験則と指針を解説します。

 

PLUMEDの設定例

MetadynamicsはPLUMEDのMETADコマンド(https://www.plumed.org/doc-v2.9/user-doc/html/_m_e_t_a_d.html?ref=molmod.dsf.unica.it)を使うことで実行できます。

以下は、PLUMEDでMetadynamicsを実行するコマンドの一例です。

# 単位設定 (長さ: Å, エネルギー: eV)
UNITS LENGTH=A ENERGY=eV

# CV定義 (原子間距離) ※原子IDは1始まり
dist: DISTANCE ATOMS=1,2

# 壁ポテンシャル (距離上限 3.5Å)
uwall: UPPER_WALLS ARG=dist AT=3.5 KAPPA=15.0 EXP=2 EPS=1 OFFSET=0

# Metadynamics設定 (300K, BIASFACTOR 80.0)
METAD ARG=dist SIGMA=0.2 HEIGHT=0.05 PACE=200 BIASFACTOR=80.0 TEMP=300.0 LABEL=metad FILE=output/HILLS

 

ポイント1: ガウス関数の高さ (HEIGHT)

積み上げるガウス型ポテンシャルの初期の高さを設定します。

  • 推奨値: 

    • 一般的には、熱揺らぎ程度の高さから検討を始めるのが安全です。

    • 例えばT=300Kの場合、kT ~ 0.6 kcal/mol ~ 0.025 eVであるため、この付近の値がよく用いられます。

  • 調整の指針:

    • ガウス関数の高さが高すぎると系が不安定になりやすいため、まずは熱揺らぎ程度の小さい値から試すことが推奨されます。しかし、状態変化や反応がなかなか進まない場合は、より高い値を設定してバイアスの蓄積を速めることを検討すると良いです。

 

ポイント2. ガウス関数の幅 (SIGMA)

積み上げるガウス型ポテンシャルの幅を設定します。

  • 推奨値:

    • 経験的には 反応座標（CV: Collective Variables）が取りうる変動範囲の 5 ~ 10% 程度の大きさで設定すると良いでしょう。

      • (例1) 原子間距離を1-5Åの範囲でメタダイナミクスをする場合、ガウス関数の幅は範囲の5%の0.2Åと設定する。

      • (例2) 配位数など0-1の範囲をとるCVでメタダイナミクスをする場合、ガウス関数の幅は範囲の5%に相当する0.05と設定する。

  • 注意点:

    • ガウス関数の幅が大きすぎる場合、 細かいエネルギー地形の特徴が潰れてしまい、正確な自由エネルギー面（FES）が得られないことがあります。

    • ガウス関数の幅が小さすぎる場合、 エネルギー地形を埋めるのに膨大な計算ステップが必要となり、探索効率が低下します。

 

ポイント3. バイアスファクター (BIASFACTOR)

BIASFACTOR を指定するとWell-Tempered Metadynamicsが実行されます。Well-tempered Metadynamicsはシミュレーション中に積み上げるガウシアンの高さを徐々に減衰させることで収束性を改善させた手法です。BIASFACTORは、バイアスの積み上げとともにガウス関数の高さの減衰の速さを制御するパラメータです。

  • γが大きいほどガウス関数の高さの減衰は遅くなるため、高い自由エネルギー障壁を乗り越えやすくなります。一方でγ→∞で通常のメタダイナミクスと同等となるため、γが大きいときはエネルギー面が振動しやすく、FESの収束判定は難しくなります。

  • γが小さいほどガウス関数の高さの減衰は早くなりますが、自由エネルギー障壁を乗り越えるまでより長いシミュレーションが必要となることがあります。
     

BIASFACTORの設定の目安

乗り越えたい最大自由エネルギー障壁ΔF_{max}と温度Tとバイアスファクターγの間には

      γ = (1 + ΔT/T) ≈ ΔF_{max} / (kT) 

の関係式が成り立ちます。従って、乗り越えたい自由エネルギー障壁の高さから適切なγの値を推測することができます。例えば、最大障壁がΔF_{max} = 20 kcal/mol、温度 T= 300K (kT = 0.6 kcal/mol) の場合、

     γ = 20/ 0.6 = 33.3 ~ 30

と計算できます。安全策をとる場合はこれよりも小さい値 (γ = 10 ~ 15)程度から試すことをお勧めします。

 

ポイント4. 壁ポテンシャル (UPPER_WALLS / LOWER_WALLS)

メタダイナミクスを効率的に実行するために、反応座標の探索範囲を限定する壁ポテンシャルを設定することが効果的なことがあります。

メタダイナミクスは、ガウス関数を積み上げることですでに探索した場所のエネルギーを底上げし、未知の領域へと系を押し出す手法です。しかし、CVの可動域に制限がない場合、例えば、結合解離の計算などで、原子・分子間距離が離れすぎると、とり得る配置の数（エントロピー）が増大し、メタダイナミクスシミュレーション中に系が結合状態に戻ってこられなくなる問題が生じます。

これらを防ぐために、UPPER_WALLS（上限）や LOWER_WALLS（下限）を設置し、系を自由エネルギーの探索範囲に適した範囲に制限することが効果的となります。

関連

関連記事