Stochastic Surface Walking (SSW) 法の理論的枠組みとポテンシャルエネルギー曲面探索における応用

分子動力学 (MD) 法: 有限温度でのサンプリングに広く用いられますが、タンパク質のフォールディングやカーボンナノチューブの成長のような高次元かつ高障壁な系では、予測能力が著しく低下します。
Simulated Annealing (SA): 温度を操作して障壁を越える手法ですが、自由エネルギー的に不利な極小点や、高温で存在しない極小点の探索には効率的でない場合があります。
Basin Hopping (BH) 法: Walesらによって体系化されたこの手法は、PESを階段状に変形（Basinの底のエネルギーにマッピング）させることでGM探索に極めて高い効率を示します [cite: 44, 45, 46]。しかし、探索過程で経路情報を破棄するため、反応経路に関する情報は得られません。
TS探索法: 固有ベクトル追跡（Eigenvector-following）などは反応経路を特定できますが、多くの場合ヘシアン（エネルギーの二階微分）の計算が必要となり、大規模系への適用には計算コストの課題が伴います。

1.2 SSW法の位置づけ#

ShangとLiu (2013) が提案した Stochastic Surface Walking (SSW) 法は、これらの課題を克服するために設計されたアルゴリズムです。SSW法は、バイアスポテンシャル駆動型ダイナミクス (Bias-Potential-Driven Dynamics) と Metropolis モンテカルロ (MC) 法 を組み合わせたものであり、以下の特長を有します。

ヘシアン不要: 勾配（力）のみを用いて構造変形を行うため、大規模系への適用が可能。
経路情報の保持: GM探索だけでなく、準安定状態間を結ぶ反応経路や障壁の情報を探索軌跡として保持する。
Unbiasedな探索: 事前の反応座標や生成物の情報を必要とせず、ランダムな初期構造から探索を開始できる。

本稿では、SSW法の数学的定式化、アルゴリズムの詳細、および検証結果について記述します。

2. 数学的枠組みとアルゴリズム#

SSW法の1ステップ（MCステップ）は、大きく分けて (i) Climbing（登山過程）、(ii) Relaxation（緩和過程）、(iii) Metropolis MC（採択判定）の3段階で構成されます。特にClimbing過程における方向生成とバイアスポテンシャルの付加が本手法の核心です。

2.1 Climbing過程：バイアスポテンシャルによる構造変形#

現在着目している極小点（Minimum）を $\mathbf{R}_i^m$ とします（ $i$ はMCステップ数）。Climbing過程では、この安定点からポテンシャルエネルギー曲面を「登り」、高エネルギー配置 $\mathbf{R}_i^H$ へと系を誘導します。

2.1.1 ランダム探索方向の生成 ( $\mathbf{N}_i^0$ )#

SSWでは、完全にランダムな方向へ変形させるのではなく、物理的に意味のある変形モードを効率的に探索するために、ソフトな大域的変位（Global move）とハードな局所的変位（Local move）を混合した初期ベクトル $\mathbf{N}_i^0$ を生成します。

\mathbf{N}_i^0 = \frac{\mathbf{N}_i^g + \lambda \mathbf{N}_i^l}{|\mathbf{N}_i^g + \lambda \mathbf{N}_i^l|} \quad (1)

ここで、 $\lambda$ は混合パラメータであり、各ステップにおいて 0.1 ～ 1.5 の範囲のランダムな値をとります。

$\mathbf{N}^g$ (Global component): 全原子に対して、300 KでのMaxwell-Boltzmann速度分布に従うランダムな正規化ベクトルを与えます。これは系全体の協奏的な動き（Soft global move）を表します。
$\mathbf{N}^l$ (Local component): 特定の化学結合の形成・解離を促すモードです。非隣接原子対（A, B、距離 > 3 Å）をランダムに選び、その結合形成方向へ変位させます。

具体的には、原子Aの座標 $\mathbf{q}_A$ と原子Bの座標 $\mathbf{q}_B$ を用いて、以下のようなベクトル成分を持ちます（正規化前）。

\mathbf{N}^l \propto (\mathbf{q}_A - \mathbf{q}_B) \quad (2)

この局所モードの導入により、共有結合の組み換えを伴う化学反応の探索効率が向上します。

2.1.2 Biased Dimer Rotationによる方向の洗練#

初期ベクトル $\mathbf{N}_i^0$ はランダム性が強すぎるため、そのままではエネルギーの高い「壁」に衝突する可能性があります。そこで、Biased Dimer Rotation 法を用いて、局所的な曲率（Hessianの固有ベクトル情報）を反映した「登りやすい方向」へベクトルを修正します。

Dimer法（Henkelman & Jónsson, 1999）は、2つのイメージ $\mathbf{R}_0, \mathbf{R}_1$ （距離 $\Delta R$ ）を用いて、ヘシアンを計算せずにその最低固有モード（ソフトモード）を見つける手法です。SSWでは、これを改良し、初期方向 $\mathbf{N}_i^0$ から大きく逸脱しない範囲でソフトモードを探すバイアス項 $V_N$ を導入しています。

V_{R1} = V_{real} + V_N \quad (5)

V_N = -\frac{a}{2} [(\mathbf{R}_1 - \mathbf{R}_0) \cdot \mathbf{N}_i^0]^2 = -\frac{a}{2} (\Delta R \cdot \mathbf{N}_t \cdot \mathbf{N}_i^0)^2 \quad (6)

ここで $a$ は拘束の強さを決めるパラメータです。このバイアスにより、初期のランダム方向の成分を保ちつつ、局所的なPESの谷に沿った方向へ修正されたベクトル $\mathbf{N}_i'$ が得られます。

2.1.3 ガウス型バイアスポテンシャルの付加#

方向 $\mathbf{N}_i'$ が定まると、その方向に沿って系を変位させます。単に変位させるのではなく、通った経路を埋めるようにガウス型のバイアスポテンシャル $v_n$ を次々と追加していきます。これにより、元の極小点に戻ることなく、高エネルギー領域（遷移状態付近）へと押し上げられます。

修正されたPES（Modified PES）は以下のように記述されます。

V_{m-to-H} = V_{real} + \sum_{n=1}^{H} v_n \quad (7)

ここで、各バイアスポテンシャル $v_n$ は以下のガウス関数です。

v_n = w_n \times \exp \left[ - \frac{((\mathbf{R}^t - \mathbf{R}_t^{n-1}) \cdot \mathbf{N}_i^n)^2}{2 \times ds^2} \right] \quad (8)

$w_n$ : ガウス関数の高さ。
$ds$ : ガウス関数の幅（変位のステップサイズ）。典型的には 0.2 ～ 0.6 Å 。
$\mathbf{R}_t^{n-1}$ : $n-1$ 番目のバイアス中心。
$H$ : ガウス関数の最大数（典型的には14程度）。

このプロセスは、Metadynamics法におけるヒストリー依存電位の付加に類似していますが、SSWでは全空間を埋めるのではなく、特定の探索経路に沿って局所的にポテンシャルを持ち上げる点が異なります。

2.2 RelaxationとMetropolis MC#

Climbing過程は、最大ステップ数 $H$ に達するか、あるいは修正PES上での局所最適化によって新たな極小点が見つかるまで続けられます。高エネルギー点 $\mathbf{R}_i^H$ に到達した後、追加した全バイアスポテンシャルを除去し、真のPES ( $V_{real}$ ) 上で構造最適化（Relaxation）を行います。これにより、新たな極小点 $\mathbf{R}_i^{mt}$ が得られます。

最後に、Metropolis基準に従って、この新構造 $\mathbf{R}_i^{mt}$ を受理するか判定します。

P = \begin{cases} \exp \left[ (E(\mathbf{R}_i^{m}) - E(\mathbf{R}_i^{mt})) / RT \right] & \text{when } E(\mathbf{R}_i^{mt}) > E(\mathbf{R}_i^{m}) \\ 1 & \text{otherwise} \end{cases} \quad (9)

ここで $E(\mathbf{R})$ は構造 $\mathbf{R}$ のエネルギー、 $R$ は気体定数（または $k_B$ ）、 $T$ は温度です。受理されれば次のステップの始点とし、拒絶されれば元の構造に戻ります。この一連の手続きにより、詳細釣り合い条件を満たしつつ、PES上を確率的に「歩き回る（Walking）」ことが可能となります。

3. 実証実験と結果#

ShangとLiuは、この手法の有効性を検証するために、小分子からクラスター、ナノ構造体に至るまで多様な系に適用しました。

3.1 ブタジエン ( $C_4H_6$ ) の配座探索#

単純な有機分子であるブタジエンの異性化反応において、SSW法は以下の挙動を示しました。

異性体の発見: 90ステップの単一走査において、最も安定な trans-ブタジエンから出発し、cis-ブタジエン、シクロブテン（環化反応）、および高エネルギーラジカル構造の計4つの異性体を発見しました。
障壁の推定: 探索軌跡上の最高エネルギー点 ( $E_{max}$ ) を用いて反応障壁を見積もったところ、cis-trans 異性化で0.36 eV、閉環反応で1.75 eVという値が得られました。これらはTS探索によって厳密に求められた値（それぞれ0.33 eV, 1.56 eV）と良い一致を示しています。ただし、SSW法ではTSを明示的に探索するわけではなく、この $E_{max}$ はあくまで障壁の近似値であることに留意が必要です。

3.2 Lennard-Jones (LJ) クラスターの大域的最適化#

$LJ_{75}$ クラスターは、PESが「二重の漏斗（Double Funnel）」構造を持つことで知られる難問です。大域的安定構造（GM）であるMarks正十面体構造と、第二安定構造（SLM）であるMackey正二十面体構造の間には、大きな構造変化と高いエネルギー障壁が存在します。

探索効率: SSW法は、SLMからスタートしてGMへ到達する遷移を効率的に発見しました。平均 $5.3 \times 10^4$ ステップで50%の確率でGMを発見し、高い障壁乗り越え能力を示しました。
経路のバイパス: 興味深いことに、SSWが見つけた経路は、高エネルギーの極小点を多数経由する従来の経路とは異なり、多くの極小点をバイパスして直接的にGMのベイスンへ遷移する傾向が見られました。これはMetropolis MCの評価回数を減らし、計算効率を向上させる要因となります。

3.3 Morseクラスターとナノヘリックス#

より相互作用距離が短く、PESが「波打った（Corrugated）」形状を持つMorseポテンシャルクラスター（ $M_{80}$ など）においても、SSW法は高いGM発見率を記録しました。これは、Biased Dimer法によるソフトモード追跡が、複雑な曲面での方向決定に有効であることを示しています。

また、Boerdijk-Coxeter-Bernal (BCB) ナノヘリックス（ $LJ_{74}$ ）の崩壊過程のシミュレーションでは、ポテンシャル障壁の上限値 ( $E_r$ ) を制御パラメータとして用いることで、準安定中間体（Intermediate）の寿命（滞在時間）を解析し、反応のカイネティクスに関する知見を引き出すことに成功しています。

4. 考察：SSW法の優位性と意義#

4.1 “Soft” なモードと “Hard” なモードの統合#

SSW法の設計思想における重要な点は、式(1)における大域的モードと局所的モードの混合、および式(6)のBiased Dimerによる修正です。従来のランダム変位のみに頼る手法では、結合長の変化などの「硬い（Hard）」モードに対してエネルギーが急上昇し、物理的に意味のある構造変化を起こす前に拒絶される傾向がありました。SSW法は、Dimer法を用いてHessianの最低固有値方向（最も曲率の小さい方向＝ソフトモード）を数値的に探索するため、高次元空間においても「谷底」に沿った、エネルギー的に無理のない変形経路を見つけ出すことができます。

4.2 反応経路の自動サンプリング#

TS探索法（NEB法やString法など）は始点と終点の構造を事前に知る必要がありますが、SSW法は一端（現在の極小点）からの探索が可能です。これにより、「予期せぬ反応経路」 や 「未知の準安定構造」 を発見する能力を持ちます。これは、触媒反応機構の解明や、新奇なアモルファス構造・結晶多形の探索において極めて強力なツールとなります。

4.3 第一原理計算との親和性#

SSW法はヘシアン（ $3N \times 3N$ 行列）の計算を必要とせず、エネルギーと力（勾配）のみで動作します。また、ステップごとの変位量 $ds$ が小さいため、電子状態計算におけるSCF（自己無撞着場）計算の収束性を損ないにくいという利点があります。このため、計算コストの高いDensity Functional Theory (DFT) を用いた第一原理計算においても、現実的な計算時間で適用可能です。

5. 今後の展望#

原著論文（2013年）の発表以降、SSW法はさらなる発展を遂げています。特に、SSW法を用いてPESを網羅的にサンプリングし、そのデータを教師データとしてニューラルネットワークポテンシャル（NNP）を構築する手法（SSW-Reaction Sampling; SSW-RSなど）が開発されています。これにより、第一原理精度の信頼性を保ちつつ、より大規模かつ長時間スケールの反応ダイナミクスシミュレーションが可能になりつつあります（注：NNPへの応用は論文発表後の発展であり、原著では将来の方向性として示唆されるにとどまっています）。

6. 結論#

ShangとLiu (2013) によって提案されたStochastic Surface Walking (SSW) 法は、バイアスポテンシャルを用いた滑らかな構造摂動と、Metropolis法による統計的サンプリングを巧みに融合させたアルゴリズムです。本手法は、事前の知識なしに複雑なPES上の大域的安定構造と反応経路を同時に探索できる「一般目的（General-purpose）」の手法として位置づけられます。特に、Biased Dimer Rotationによる方向制御とガウス関数による局所的な障壁乗り越え機構は、高次元かつ荒いPESを持つ系に対して高い頑健性を示しました。

参考文献#

[1] Shang, C.; Liu, Z.-P. “Stochastic Surface Walking Method for Structure Prediction and Pathway Searching.” J. Chem. Theory Comput. 2013, 9, 1838–1845. DOI: 10.1021/ct301010b
[2] Wales, D. J.; Doye, J. P. K. “Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms.” J. Phys. Chem. A 1997, 101, 5111.
[3] Henkelman, G.; Jónsson, H. “A Dimer Method for Finding Saddle Points on High Dimensional Potential Surfaces Using Only First Derivatives.” J. Chem. Phys. 1999, 111, 7010.
[4] Laio, A.; Parrinello, M. “Escaping Free-Energy Minima.” Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562.