Doubly Nudged Elastic Band (DNEB) 法の理論的枠組みと準ニュートン法による遷移状態探索の効率化

最終更新：2026-01-02

注意: この記事はGemini 3.0によって自動生成されたものです。内容はSemen A. TrygubenkoおよびDavid J. Walesによる原著論文 “A doubly nudged elastic band method for finding transition states” (J. Chem. Phys. 120, 2082, 2004) に基づきます。正確な学術情報については、必ず原著論文を参照してください。

序論：ポテンシャルエネルギー曲面における遷移状態探索の課題#

化学反応動力学や材料科学における原子拡散過程の理解において、ポテンシャルエネルギー曲面（Potential Energy Surface; PES）上の停留点、特に極小点（Minima）と遷移状態（Transition State; TS）を特定することは不可欠な手続きである。MurrellとLaidlerの幾何学的定義によれば、遷移状態はヘッセ行列（Hessian）がただ一つの負の固有値を持つ停留点、すなわち一次の鞍点（Saddle Point）として定義される。

遷移状態の探索手法は、単一の初期構造から探索を開始するSingle-ended法（固有ベクトル追跡法など）と、反応物と生成物の二つの極小点間を結ぶ経路を探索するDouble-ended法に大別される。Double-ended法の代表格であるNudged Elastic Band (NEB) 法は、直感的な概念と実装の容易さから広く普及している。しかし、2000年代初頭において、NEB法は強力な準ニュートン法（L-BFGS等）と組み合わせた際に数値的な不安定性を示すことが知られており、最適化効率の向上が課題となっていた。

本稿では、TrygubenkoとWales (2004) によって提案されたDoubly Nudged Elastic Band (DNEB) 法について詳述する。本手法は、標準的なNEB法の勾配定義を修正することで、L-BFGS法を用いた際の安定性を回復させ、計算効率を最大で2桁向上させるものである。

1. 歴史的背景とChain-of-States法の進化#

NEB法の位置づけを明確にするため、Double-ended法の歴史的変遷を概観する。

1.1 初期の手法：LSTとQST#

最も初期の試みとして、線形同期通過（Linear Synchronous Transit; LST）および二次同期通過（Quadratic Synchronous Transit; QST）が挙げられる。

LST: 始点と終点を結ぶ直線上で最高エネルギーを持つ構造を探索する。
QST: 反応経路を放物線で近似する。

これらの補間ベースの手法は、計算コストは低いものの、複雑な反応経路においては実際の遷移状態から大きく外れた構造を与える傾向があり、より洗練された手法の初期推定としての利用に留まることが多い。

1.2 Chain-of-Statesアプローチ#

経路を離散的な「イメージ（Image）」の連なりとして表現し、それらを同時に最適化するChain-of-States (CS) 法が登場した。初期のCS法では、イメージ間の距離を拘束条件として課したり、ペナルティ関数を用いたりしていたが、経路が曲率を持つ領域で「Corner-cutting（近道）」の問題が発生しやすかった。

1.3 Nudged Elastic Band (NEB) 法の登場#

Jónssonら（1998）によって導入されたNEB法は、イメージ間に仮想的なバネを導入しつつ、力の射影（Nudging）を行うことでCS法の問題を解決しようと試みた。 NEB法の核心は、ポテンシャル力（真の力）とバネ力の干渉を防ぐために、以下のように力の成分を分解・除去することにある。

ポテンシャル力: 経路の接線方向成分を除去する（Sliding-downの防止）。
バネ力: 経路の接線に対する垂直成分を除去する（Corner-cuttingの防止）。

しかし、この「Nudging」操作により、NEBの力場は非保存力場（Non-conservative force field）となり、対応する大域的な目的関数（Global Objective Function）が存在しなくなる。これが、ヘッセ行列の対称性を前提とする準ニュートン法との相性を悪化させる根本原因であった。

2. NEB法の数理的定式化と不安定性の起源#

ここで、標準的なNEB法の数理的構造と、Trygubenkoらが指摘した問題点を整理する。

2.1 力の分解と射影#

$N$ 個のイメージからなるバンドを考える。イメージ $i$ の座標ベクトルを $\mathbf{X}_i$ とし、ポテンシャルエネルギーを $V(\mathbf{X})$ とする。イメージ $i$ における真の勾配（ポテンシャルの勾配）を $\mathbf{g}_i = \nabla V(\mathbf{X}_i)$ 、バネによる勾配を $\tilde{\mathbf{g}}_i$ とする。経路の接線単位ベクトルを $\hat{\boldsymbol{\tau}}_i$ と定義するとき、勾配の平行成分（ $\parallel$ ）と垂直成分（ $\perp$ ）は以下のように表される。

\mathbf{g}_i^{\parallel} = (\mathbf{g}_i \cdot \hat{\boldsymbol{\tau}}_i) \hat{\boldsymbol{\tau}}_i, \quad \mathbf{g}_i^{\perp} = \mathbf{g}_i - \mathbf{g}_i^{\parallel}

標準的なNEB法における有効勾配 $\mathbf{g}^{NEB}_i$ は、以下のように定義される。

\mathbf{g}^{NEB}_i = \mathbf{g}_i^{\perp} + \tilde{\mathbf{g}}_i^{\parallel}

ここで、 $\tilde{\mathbf{g}}_i^{\parallel}$ は隣接イメージ間の距離に基づくバネ力の平行成分である。この定義により、イメージはポテンシャルの谷底を経由しつつ（ $\mathbf{g}_i^{\perp}$ の寄与）、経路上で等間隔に配置される（ $\tilde{\mathbf{g}}_i^{\parallel}$ の寄与）ことが期待される。

2.2 接線ベクトルの推定#

接線 $\hat{\boldsymbol{\tau}}_i$ の推定方法は収束性に大きく影響する。初期のNEBでは隣接する両側のイメージ（ $i+1$ と $i-1$ ）を用いていたが、これは経路が急激に曲がる箇所で「Kinks（折れ曲がり）」を生じさせることが判明した。これを改善するため、HenkelmanとJónsson（2000）は、エネルギーが高い方の隣接イメージを用いて接線を定義する手法（Upwind tangent scheme）を導入した。

\hat{\boldsymbol{\tau}}_i = \frac{\mathbf{X}_j - \mathbf{X}_i}{|\mathbf{X}_j - \mathbf{X}_i|} \quad \text{where } E_j > E_i

2.3 L-BFGS適用時の不安定性#

TrygubenkoとWalesは、標準的なNEB勾配を用いてL-BFGS（Limited-memory Broyden-Fletcher-Goldfarb-Shanno）法による最適化を行った場合、計算が不安定になり、収束しないか、あるいは誤った経路へ発散することを示した。

この不安定性の主要因は、バネ力の垂直成分 $\tilde{\mathbf{g}}_i^{\perp}$ を完全に除去している点にある。 L-BFGSなどの準ニュートン法は、勾配の履歴からヘッセ行列（曲率）を近似する。標準NEBでは、イメージが経路から外れた際に、それを元の経路へ引き戻す復元力（バネの垂直成分）が存在しないか、あるいはポテンシャルの垂直成分のみに依存することになる。経路が大きく曲がっている領域では、ポテンシャルの曲率情報と、幾何学的な拘束（バネ）の情報の不整合が生じ、近似ヘッセ行列が不適切な探索方向を導き出す結果となる。特に、L-BFGSが大きなステップをとろうとした際に、一時的な不連続性やKinksが発生し、アルゴリズムが破綻する。

3. Doubly Nudged Elastic Band (DNEB) 法#

上述の不安定性を解消するために提案されたのが、Doubly Nudged Elastic Band (DNEB) 法である。

3.1 DNEBの基本概念#

DNEBのアイデアは、標準NEBで除去されていた「バネ力の垂直成分 $\tilde{\mathbf{g}}_i^{\perp}$ 」の一部を、修正を加えた上で復元することである。これにより、イメージ列の「弾性バンド」としての性質（＝滑らかさを保とうとする性質）を部分的に復活させ、L-BFGSに対して安定な曲率情報を提供する。

単純に $\tilde{\mathbf{g}}_i^{\perp}$ を全て加えてしまうと、本来のNEBが解決したCorner-cutting問題が再発する。そこでTrygubenkoらは、ポテンシャルの勾配に影響を与えない成分のみを復元する射影を考案した。

3.2 数学的導出#

DNEBにおける勾配 $\mathbf{g}^{DNEB}_i$ は以下のように定義される。

\mathbf{g}^{DNEB}_i = \mathbf{g}_i^{\perp} + \tilde{\mathbf{g}}_i^{\parallel} + \tilde{\mathbf{g}}_i^{*}

ここで、 $\tilde{\mathbf{g}}_i^{*}$ がDNEB特有の修正項（Double Nudging項）である。この項は、バネの垂直成分 $\tilde{\mathbf{g}}_i^{\perp}$ から、さらに「ポテンシャルの垂直勾配 $\mathbf{g}_i^{\perp}$ 」方向への射影成分を除去したものとして定義される。

\tilde{\mathbf{g}}_i^{*} = \tilde{\mathbf{g}}_i^{\perp} - (\tilde{\mathbf{g}}_i^{\perp} \cdot \hat{\mathbf{g}}_i^{\perp}) \hat{\mathbf{g}}_i^{\perp}

ここで、 $\hat{\mathbf{g}}_i^{\perp}$ はポテンシャルの垂直勾配方向の単位ベクトルである（ $\mathbf{g}_i^{\perp} = 0$ の場合は $\tilde{\mathbf{g}}_i^{*} = 0$ とする）。

3.3 物理的解釈#

この射影操作の意味するところは以下の通りである。

$\mathbf{g}_i^{\perp}$ 成分の除去: バネ力が、ポテンシャルの壁を登る方向（あるいは下る方向）に寄与することを防ぐ。これにより、イメージがポテンシャルの谷底から引きずり出されることを回避し、Corner-cuttingを最小限に抑える。
残りの垂直成分の保持: ポテンシャル勾配と直交する方向（谷の断面内での横方向の揺らぎなど）に対しては、バネの復元力を作用させる。これがバンドを「ピンと張った」状態に保つ役割を果たし、L-BFGSにおける数値的な安定性に寄与する。

この「二重のナッジング（Doubly Nudged）」により、DNEBはNEBの利点（正確な経路追跡）と、従来のChain-of-States法の利点（最適化の安定性）を両立させている。

4. 最適化アルゴリズムの実装と評価#

TrygubenkoとWalesは、DNEB法の性能を評価するために、最適化アルゴリズムの詳細な検討を行っている。ここでは、比較対象となった手法と、その実装上の要点を解説する。

4.1 Quenched Velocity Verlet (QVV) とその改良#

NEBの最適化には従来、分子動力学法（MD）に基づくQuenched Velocity Verlet (QVV) 法が用いられてきた。これは、速度ベクトルが力と逆方向（エネルギーが増加する方向）を向いた瞬間に速度をゼロにする（Quenching）ことで、エネルギー散逸を模倣し、極小点へ収束させる手法である。

論文では、速度のリセット（Quenching）を行うタイミングについて詳細な比較検討が行われている。

標準的なQVV: 速度更新の都度、内積 $\mathbf{v} \cdot \mathbf{g}$ をチェックする。
Slow-response QVV (SQVV): 座標更新の直後にQuenching判定を行う。具体的には、座標更新後の新しい勾配ではなく、1ステップ前の速度情報を用いて判定を行う形となる。

Trygubenkoらは、SQVV（座標更新後のQuenching）が最も安定性が高く、バネ定数 $k_{spr}$ の値に対してロバストであることを示した。従来のQVVは、バネ定数が大きい場合に振動的な挙動を示しやすく、収束が遅くなる傾向があった。

4.2 Limited-memory BFGS (L-BFGS) の適用#

L-BFGSは、ヘッセ行列全体を保持せず、直近 $m$ ステップの変位と勾配変化の履歴から逆ヘッセ行列との積を計算する。大規模系に適しているが、上述の通り標準NEBでは不安定であった。

DNEB法においてL-BFGSを適用する際、以下の実装上の工夫がなされた。

Global L-BFGSアプローチ: 全イメージの全座標を一つの巨大なベクトル（次元数 $N_{images} \times 3N_{atoms}$ ）として扱い、L-BFGSを適用する。
ステップサイズの制限: 信頼領域法（Trust Region Method）の考え方に従い、各イメージの最大移動量を制限する。これにより、バンドの絡まりや極端な変形を防ぐ。
初期ヘッセ行列のスケーリング: 単位行列ではなく、適切なスカラー倍を用いることで初期の探索効率を改善する。

4.3 性能比較：DNEB/L-BFGS vs NEB/SQVV#

Lennard-Jonesクラスター（ $LJ_7, LJ_{38}, LJ_{75}$ ）の異性化反応をテストケースとして、両手法の性能比較が行われた。

収束速度: DNEB/L-BFGSは、最適化完了までの反復回数（および勾配評価回数）において、NEB/SQVVと比較して1桁から2桁（約10倍〜100倍）の高速化を達成した。
ロバスト性: 標準NEBではL-BFGSが発散するような複雑な経路においても、DNEBは安定して収束した。
バネ定数依存性: NEB/SQVVはバネ定数の設定に敏感で、不適切な値では収束が著しく悪化したが、DNEB/L-BFGSは広い範囲のバネ定数に対して安定した性能を示した。

5. 高度な接続アルゴリズムと実用例#

DNEB法の高速性は、単一の遷移状態探索だけでなく、複雑な多段階反応経路の全貌解明（Discrete Path Sampling; DPS）においても威力を発揮する。

5.1 逐次的な経路構築戦略#

$LJ_{38}$ や $LJ_{75}$ のような二重漏斗（Double-funnel）型エネルギーランドスケープを持つ系では、反応物と生成物の構造が大きく異なり、その間には多数の中間体と遷移状態が存在する。これらを一度のNEB計算で繋ぐことは困難である。

Trygubenkoらは、DNEBを用いて以下の手順で経路を逐次的に構築するアルゴリズムを提案した。

初期探索: 始点と終点の間でDNEBを実行し、最もエネルギーの高いイメージを遷移状態候補として抽出する。
TSの精密化: 抽出された候補構造に対し、固有ベクトル追跡法（Eigenvector-following）を適用して正確な遷移状態へ収束させる。
IRC計算: 得られたTSから最急降下経路（Intrinsic Reaction Coordinate; IRC）を計算し、新たな中間安定点を見つける。
データベース更新: 未接続の極小点ペアに対して、距離が近い順に再度DNEBを適用し、経路網（Transition Network）を拡張していく。

5.2 結果の実証#

この自動化された手順により、 $LJ_{38}$ クラスターにおける秩序-無秩序転移や、 $LJ_{75}$ における複雑な協同的再配置過程の全容が解明された。DNEB/L-BFGSの高い効率性が、数千〜数万回に及ぶ遷移状態探索を現実的な計算時間で実行することを可能にしたのである。

6. 結論と展望#

TrygubenkoとWales (2004) によるDNEB法の提案は、Chain-of-States法における最適化アルゴリズムの問題に、幾何学的な視点から一つの解を与えた。

数学的貢献: バネ力の射影方法を再考し、ポテンシャル勾配と干渉しない範囲でバネの復元力を維持する「二重ナッジング（Double Nudging）」の定式化（Eq. 13）を行った。
実利的成果: これにより、強力な大域的最適化手法であるL-BFGS法の利用が可能となり、従来のMDベースの手法（SQVV）と比較して劇的な計算コストの削減を実現した。
歴史的意義: 本研究は、その後のSheppardら（2008）によるGlobal L-BFGSの再評価や、現在の主要な第一原理計算コード（VASP, VTST Tools等）におけるNEB実装の基礎となる重要な知見を提供した。

DNEB法は、現代の計算化学においても、複雑な反応経路を高速かつロバストに探索するための標準的なツールセットの一つとして位置づけられている。特に、ポテンシャル曲面が平坦で探索が困難な系や、多数の自由度を持つ凝縮系において、その有効性は揺るぎないものである。

付録：DNEB勾配計算のPython擬似コード#

以下に、DNEB法における勾配計算部分の概念的な実装をPython (NumPy) 形式で示す。

import numpy as np

def calculate_dneb_forces(images, energies, potential_forces, spring_constant):
    """
    DNEB法に基づく力の計算 (概念実装)
    
    Parameters:
    -----------
    images : np.ndarray (N_images, N_atoms, 3)
        各イメージの座標
    energies : np.ndarray (N_images,)
        各イメージのポテンシャルエネルギー
    potential_forces : np.ndarray (N_images, N_atoms, 3)
        各イメージの物理的な力 (-grad V)
    spring_constant : float
        バネ定数 k
        
    Returns:
    --------
    dneb_forces : np.ndarray (N_images, N_atoms, 3)
        DNEB法によって補正された力
    """
    n_images = len(images)
    dneb_forces = np.zeros_like(images)
    
    # 接線ベクトルの計算 (Upwind scheme)
    tangents = np.zeros_like(images)
    for i in range(1, n_images - 1):
        v1 = images[i] - images[i-1]
        v2 = images[i+1] - images[i]
        e0 = energies[i-1]
        e1 = energies[i]
        e2 = energies[i+1]
        
        # エネルギー形状に応じた接線の切り替え
        if e1 > e0 and e1 > e2: # 極大点付近
             tangents[i] = v2 if abs(e2-e1) > abs(e0-e1) else v1
        elif e1 < e0 and e1 < e2: # 極小点付近
             tangents[i] = v2 if abs(e2-e1) > abs(e0-e1) else v1
        else:
            if e2 > e0:
                tangents[i] = v2
            else:
                tangents[i] = v1
        
        # 正規化
        norm = np.linalg.norm(tangents[i])
        if norm > 0:
            tangents[i] /= norm

    # 力の射影とDNEB項の計算
    for i in range(1, n_images - 1):
        # 1. バネ力 (Spring Force) の計算
        # ||R_{i+1} - R_i|| - ||R_i - R_{i-1}||
        dist_next = np.linalg.norm(images[i+1] - images[i])
        dist_prev = np.linalg.norm(images[i] - images[i-1])
        f_spring_scalar = spring_constant * (dist_next - dist_prev)
        f_spring_parallel = f_spring_scalar * tangents[i]
        
        # バネ力のベクトルそのもの（隣接間距離ベクトルに基づく）
        # 簡単のため接線方向成分が支配的と仮定して近似することもあるが、
        # DNEBでは垂直成分が必要なため、厳密には以下のようなベクトル差分を考慮する
        f_spring_vec = spring_constant * ((images[i+1] - images[i]) + (images[i-1] - images[i]))
        
        # 2. ポテンシャル力 (Potential Force) の分解
        f_pot = potential_forces[i]
        f_pot_parallel_mag = np.vdot(f_pot, tangents[i])
        f_pot_parallel = f_pot_parallel_mag * tangents[i]
        f_pot_perpendicular = f_pot - f_pot_parallel # g_perp
        
        # ポテンシャル垂直方向の単位ベクトル g_hat_perp
        g_perp_norm = np.linalg.norm(f_pot_perpendicular)
        if g_perp_norm > 1e-10:
            g_hat_perp = f_pot_perpendicular / g_perp_norm
        else:
            g_hat_perp = np.zeros_like(f_pot_perpendicular)

        # 3. バネ力の垂直成分 (g_tilde_perp)
        # 注意: 論文中の記法では勾配 g = -F なので符号に注意。
        # ここでは力 F を扱っている。
        f_spring_perpendicular = f_spring_vec - np.vdot(f_spring_vec, tangents[i]) * tangents[i]

        # 4. DNEB補正項 (Eq. 13 in Trygubenko2004)
        # tilde_g* = tilde_g_perp - (tilde_g_perp . g_hat_perp) * g_hat_perp
        # 力ベースの記法に直すと:
        # F* = F_spring_perp - (F_spring_perp . (-g_hat_perp)) * (-g_hat_perp)
        # 内積部分の符号は相殺されるため形は同じ
        projection = np.vdot(f_spring_perpendicular, g_hat_perp)
        f_dneb_correction = f_spring_perpendicular - projection * g_hat_perp

        # 5. 最終的なDNEB力
        # F_DNEB = F_pot_perp + F_spring_parallel + F_DNEB_correction
        dneb_forces[i] = f_pot_perpendicular + f_spring_parallel + f_dneb_correction

    return dneb_forces

参考文献#

S. A. Trygubenko and D. J. Wales, A doubly nudged elastic band method for finding transition states, J. Chem. Phys. 120, 2082 (2004).
G. Henkelman and H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113, 9978 (2000).
H. Jónsson, G. Mills, and K. W. Jacobsen, Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions, in Classical and Quantum Dynamics in Condensed Phase Simulations, World Scientific (1998).
D. Sheppard, R. Terrell, and G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys. 128, 134106 (2008).