Nudged Elastic Band (NEB) 法における大域的最適化アルゴリズムの包括的評価：Global L-BFGS法の数理的構造と実装論

最終更新：2026-01-02

注意: この記事はGemini 3.0によって自動生成されたものです。正確な学術情報については、必ず末尾に記載されたSheppard, Henkelmanらの原著論文（J. Chem. Phys., 2008）および関連する専門書をご確認ください。

序論：反応経路探索における計算効率の課題#

化学反応、触媒作用、および材料中の拡散現象といった原子レベルのダイナミクスを理解する上で、ポテンシャルエネルギー曲面（Potential Energy Surface; PES）上の**最小エネルギー経路（Minimum Energy Path; MEP）**を特定することは、現代計算化学における中心的課題の一つである。MEP上の最高点である鞍点（Saddle Point）は遷移状態に対応し、そのエネルギー障壁は反応速度定数を支配する主要因となる。

MEPを探索するための手法として、Nudged Elastic Band (NEB) 法（Jónsson et al., 1998; Henkelman et al., 2000）は、その堅牢性と直感的な概念から広く採用されている。NEB法は、反応物（Initial State）と生成物（Final State）を結ぶ経路を、複数の離散的な「イメージ（Image）」の連なりとして表現し、イメージ間を仮想的なバネで連結することで経路全体の最適化を行うChain-of-States法の一種である。

しかし、NEB計算、特に第一原理計算（DFT）を用いた高精度なシミュレーションにおいては、計算コストが大きな障壁となる。各イメージに働く力（NEB力）をゼロに収束させるプロセスは、高次元空間における非線形最適化問題に帰着される。従来、この最適化には最急降下法（Steepest Descent）や、分子動力学（MD）に基づいたQuick-min法などが標準的に用いられてきた。これらは実装が容易で安定している反面、収束速度は一次（Linear）に留まり、ポテンシャル曲面が複雑な形状を持つ場合や、高い精度が求められる場合には、膨大な回数のエネルギー・勾配評価を必要とする傾向があった。

2008年、Sheppard、Terrell、Henkelmanは、構造最適化の分野で実績のある準ニュートン法、特にL-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) 法をNEB法に適用し、その性能を体系的に評価した。彼らの研究は、NEB力が保存力場ではない（対応するポテンシャル関数が存在しない）という理論的な懸念にもかかわらず、系全体を一つの巨大なベクトルとして扱う「Global L-BFGS」アプローチが、従来手法と比較して優れた収束性能を示すことを実証した。

本稿では、NEB力の数理的特性、各種最適化アルゴリズムの理論的背景、そしてGlobal L-BFGS法のアルゴリズム構造とその動作原理について、数理的な導出と具体的な言語化を通じて詳細に解説する。

1. NEB法の数理的枠組みと力の特性#

最適化アルゴリズムの議論に入る前に、NEB法における「力」の定義とその特異な性質について整理する。

1.1 NEB力の定義#

NEB法において、経路を構成する $N$ 個の中間イメージ $\mathbf{R}_1, \dots, \mathbf{R}_N$ に作用する力 $\mathbf{F}_i^{NEB}$ は、ポテンシャルエネルギー曲面からの物理的な力と、イメージ間隔を制御する人工的なバネ力を、経路の接線方向に基づいて射影・合成することで定義される。

\mathbf{F}_i^{NEB} = \mathbf{F}_i^{\perp} + \mathbf{F}_i^{s, \parallel}

ここで、垂直成分 $\mathbf{F}_i^{\perp}$ は、ポテンシャル力 $-\nabla V(\mathbf{R}_i)$ のうち、経路の接線ベクトル $\hat{\boldsymbol{\tau}}_i$ に垂直な成分である。これはイメージをMEP上へと導く役割を果たす。

\mathbf{F}_i^{\perp} = -\nabla V(\mathbf{R}_i) + (\nabla V(\mathbf{R}_i) \cdot \hat{\boldsymbol{\tau}}_i) \hat{\boldsymbol{\tau}}_i

一方、平行成分 $\mathbf{F}_i^{s, \parallel}$ は、隣接イメージとの間のバネ力の接線方向成分のみを取り出したものである。これはイメージを経路に沿って等間隔（あるいは特定の密度分布）に配置する役割を果たす。

\mathbf{F}_i^{s, \parallel} = k (|\mathbf{R}_{i+1} - \mathbf{R}_i| - |\mathbf{R}_i - \mathbf{R}_{i-1}|) \hat{\boldsymbol{\tau}}_i

この「Nudging（ナッジング）」と呼ばれる射影操作により、バネ力が経路を直線化してしまう「Corner-cutting」問題と、ポテンシャル力がイメージを安定点へ引きずり落とす「Sliding-down」問題を同時に回避している。

1.2 非保存力場としての課題#

通常の構造最適化において、力 $\mathbf{F}$ はポテンシャルエネルギー $E$ の負の勾配 $\mathbf{F} = -\nabla E$ として与えられる。このような場は保存力場と呼ばれ、力の微分行列（ヘッセ行列またはヤコビアン）は対称行列となる（ $\partial F_i / \partial R_j = \partial F_j / \partial R_i$ ）。

しかし、NEB力 $\mathbf{F}^{NEB}$ は、垂直成分と平行成分を人工的に組み合わせたものであり、一般にこれを導くスカラーポテンシャル関数 $\mathcal{V}_{eff}$ は存在しない。すなわち、NEB力場は**非保存力場（Non-conservative force field）**である。

\frac{\partial F_i^{NEB}}{\partial R_j} \neq \frac{\partial F_j^{NEB}}{\partial R_i}

多くの高度な最適化アルゴリズム、特にBFGS法などの準ニュートン法は、ヘッセ行列の対称性を前提として逆ヘッセ行列を更新していく。したがって、非保存力場であるNEB力に対してこれらの手法を適用することは、数学的な厳密性の観点からは自明ではない。Sheppardらの研究の核心は、この理論的な不整合が実用上の収束性にどのような影響を与えるか（あるいは与えないか）を検証した点にある。

2. Global L-BFGS法のアルゴリズムと言語化#

本論文で最も推奨されている手法が Global L-BFGS である。これは、NEB計算に関与する全イメージの座標変数を一つの巨大な最適化問題として扱い、L-BFGS法を適用するアプローチである。

2.1 アルゴリズムの概念的構造#

L-BFGS（Limited-memory BFGS）は、ヘッセ行列の逆行列を直接保持するのではなく、直近の数ステップ分の「変位」と「力の変化」の履歴を用いて、逆ヘッセ行列とベクトルとの積を近似的に計算する手法である。Global L-BFGSでは、この手法をNEB系全体に適用する。

具体的な処理手順を以下に示す。

ステップ1: 大域的ベクトルの構築 (Flattening)#

NEB計算における全イメージ（通常は両端の固定点を除く $N$ 個の中間イメージ）の座標を、一本の巨大なベクトル $\mathbf{X}$ に結合する。

各イメージ $\mathbf{R}_i$ は $3 \times N_{atoms}$ 次元のベクトルである。
これらを連結し、次元数 $D = 3 \times N_{atoms} \times N_{images}$ の大域座標ベクトル $\mathbf{X}$ を作成する。

同様に、各イメージで計算されたNEB力 $\mathbf{F}_i^{NEB}$ も連結し、同次元の大域力ベクトル $\mathbf{F}_{total}$ を作成する。

ステップ2: 履歴情報の更新#

現在のステップ $k$ における座標 $\mathbf{X}_k$ と力 $\mathbf{F}_{total, k}$ 、および前回のステップ $k-1$ の情報を用いて、以下の二つの差分ベクトルを計算し、履歴リスト（メモリ）に保存する。

変位ベクトル $\mathbf{s}_{k-1} = \mathbf{X}_k - \mathbf{X}_{k-1}$
力の変化ベクトル $\mathbf{y}_{k-1} = (-\mathbf{F}_{total, k}) - (-\mathbf{F}_{total, k-1})$ $y_{k - 1} = (- F_{t o t a l, k}) - (- F_{t o t a l, k - 1})$
- ここで、力は負の勾配として扱われるため、符号に注意が必要である。

メモリには直近 $m$ 個（通常 $m=5 \sim 20$ ）のペア $(\mathbf{s}, \mathbf{y})$ のみを保持し、古いものは破棄する。

ステップ3: 探索方向の決定 (Two-loop Recursion)#

L-BFGSの核心部分である。現在の勾配（負の力 $-\mathbf{F}_{total, k}$ ）に対し、近似的な逆ヘッセ行列 $H_k$ を掛け合わせたベクトル、すなわちニュートンステップ方向 $\mathbf{d}_k = -H_k (-\mathbf{F}_{total, k})$ を計算する。巨大行列 $H_k$ を明示的に構成せず、以下の「Two-loop recursion」アルゴリズムを用いて計算する。

Backward Loop: 現在の勾配ベクトル $\mathbf{q} = -\mathbf{F}_{total, k}$ を初期値とし、最新の履歴から過去に向かって順に、ベクトル $\mathbf{s}_i, \mathbf{y}_i$ との内積を用いて係数 $\alpha_i$ を計算し、 $\mathbf{q}$ を更新していく。
- $\alpha_i = \rho_i \mathbf{s}_i^T \mathbf{q}$
- $\mathbf{q} \leftarrow \mathbf{q} - \alpha_i \mathbf{y}_i$
- ここで $\rho_i = 1 / (\mathbf{y}_i^T \mathbf{s}_i)$
Scaling: 中心となる初期行列 $H_0$ を、直近の曲率情報を用いてスケーリングする。
- $\mathbf{r} = H_0 \mathbf{q} = \frac{\mathbf{s}_{k-1}^T \mathbf{y}_{k-1}}{\mathbf{y}_{k-1}^T \mathbf{y}_{k-1}} \mathbf{q}$
Forward Loop: 最も古い履歴から最新のものに向かって順に、 $\mathbf{r}$ を更新していく。
- $\beta = \rho_i \mathbf{y}_i^T \mathbf{r}$
- $\mathbf{r} \leftarrow \mathbf{r} + \mathbf{s}_i (\alpha_i - \beta)$

最終的に得られた $\mathbf{r}$ が、求める探索方向 $\mathbf{d}_k$ となる。

ステップ4: 直線探索 (Line Search)#

決定された方向 $\mathbf{d}_k$ に沿って、どれだけ進むべきか（ステップサイズ $\lambda$ ）を決定する。NEB力はポテンシャルの勾配ではないため、エネルギーの減少を基準にする通常の直線探索は厳密には適用できないが、力の大きさ（二乗ノルム）の減少や、射影された力の方向との整合性を基準に行う。多くの場合、バックトラッキング法などが用いられる。

ステップ5: 座標の更新と再配置 (Reshaping)#

更新された大域座標ベクトル $\mathbf{X}_{k+1} = \mathbf{X}_k + \lambda \mathbf{d}_k$ を、元のイメージごとの座標 $\mathbf{R}_1, \dots, \mathbf{R}_N$ に分解（Reshape）し、次のステップの力計算へ進む。

2.2 Global L-BFGSの実装コード例#

以下に、PythonとNumPyを用いたGlobal L-BFGS NEBの概念的な実装を示す。これは実際の計算コードの骨格となるものである。

import numpy as np

class GlobalLBFGS_NEB:
    def __init__(self, initial_images, calc_force_func, memory_size=10):
        """
        initial_images: (N_images, N_atoms, 3) の形状を持つ初期座標配列
        calc_force_func: 現在のイメージを受け取り、NEB力を返す関数
        memory_size: L-BFGSの履歴保持数 (m)
        """
        self.images = initial_images
        self.calc_force = calc_force_func
        self.m = memory_size
        self.s_list = [] # 変位ベクトルの履歴
        self.y_list = [] # 力の変化ベクトルの履歴
        self.rho_list = [] # スケーリング係数の履歴

    def flatten(self, images):
        """(N_images, N_atoms, 3) -> (D,) の1次元ベクトルへ変換"""
        return images.flatten()

    def reshape(self, flat_vec):
        """(D,) -> (N_images, N_atoms, 3) の形状へ戻す"""
        return flat_vec.reshape(self.images.shape)

    def step(self):
        """最適化の1ステップを実行"""
        # 1. 現在のNEB力を計算
        forces = self.calc_force(self.images)
        
        # 2. 座標と力をFlatten
        x_k = self.flatten(self.images)
        f_k = self.flatten(forces) # NEB力
        grad_k = -f_k              # 勾配として扱う

        # 3. 探索方向の決定 (L-BFGS Two-loop recursion)
        d_k = self.get_lbfgs_direction(grad_k)

        # 4. ステップサイズの決定 (簡易的な固定ステップあるいは直線探索)
        # ここでは単純化のため固定ステップとするが、実際はLine Searchが必要
        alpha = 0.1 
        x_new = x_k + alpha * d_k

        # 5. 履歴の更新
        if len(self.s_list) > 0:
            # 前回の力などはクラスメンバとして保持しておく必要がある
            # ここでは概念のみ示す
            pass 

        # 実際の更新処理 (s_k, y_kの計算と保存)
        # s_k = x_new - x_old
        # y_k = grad_new - grad_old
        # if s_k @ y_k > 0:  # 曲率条件のチェック
        #     self.s_list.append(s_k)
        #     self.y_list.append(y_k)
        #     ...

        # 6. 座標更新
        self.images = self.reshape(x_new)
        return np.linalg.norm(f_k) # 収束判定用の残差

    def get_lbfgs_direction(self, q):
        """Two-loop recursionの実装"""
        if not self.s_list:
            return -q # 履歴がない場合は最急降下方向

        q = q.copy()
        alpha_list = []
        
        # Backward loop
        for i in range(len(self.s_list) - 1, -1, -1):
            s = self.s_list[i]
            y = self.y_list[i]
            rho = self.rho_list[i]
            alpha = rho * np.dot(s, q)
            alpha_list.append(alpha)
            q -= alpha * y

        # Scaling (H_0)
        s_last = self.s_list[-1]
        y_last = self.y_list[-1]
        gamma = np.dot(s_last, y_last) / np.dot(y_last, y_last)
        r = gamma * q

        # Forward loop
        alpha_list.reverse() # 順序を戻す
        for i in range(len(self.s_list)):
            s = self.s_list[i]
            y = self.y_list[i]
            rho = self.rho_list[i]
            alpha = alpha_list[i]
            beta = rho * np.dot(y, r)
            r += s * (alpha - beta)

        return -r # 探索方向 (d_k = -H * grad)

この実装における最大のポイントは、flatten と reshape によって、NEBの全イメージを一つのベクトルとして抽象化している点である。これにより、L-BFGSのロジック自体は通常の構造最適化と全く同じものを使用でき、かつイメージ間の相互作用（バネ力）がヘッセ行列の非対角成分として暗黙的に学習されることになる。

3. なぜGlobal L-BFGSは有効なのか？：物理的・数理的考察#

非保存力場に対するBFGS法の適用は、数学的には保証されないにもかかわらず、Sheppardらの結果は圧倒的な成功を示している。その理由について、数理的背景と物理的直観の両面から深く掘り下げる。

3.1 「硬いモード」と「柔らかいモード」の分離#

NEB法には、ポテンシャル曲面由来の物理的な力と、イメージ間隔を保つための人工的なバネ力という、性質の異なる二つの力が混在している。

バネ定数 $k$ の影響: 経路の解像度を上げるためにバネ定数を大きくすると、系には非常に高い周波数の振動モード（バネの伸縮）が生じる。これを「硬い（Stiff）」系と呼ぶ。

最急降下法やMD法の限界: これらの手法は、系の最も高い振動数（硬いモード）によってステップ幅（タイムステップ）が制限される。バネが硬いと、ポテンシャル曲面上の移動（柔らかいモード）に必要な大きなステップをとることができず、振動しながら極めてゆっくりとしか進まない。

Global L-BFGSの役割 (Preconditioning)における逆ヘッセ行列 $H$ の更新は、変数のスケーリングを自動的に調整する効果（Preconditioning）を持つ。アルゴリズムは履歴情報から、「バネ方向は力が急激に変化する（曲率が大きい＝硬い）」一方で、「経路の接線に垂直な方向は力が緩やかに変化する（曲率が小さい＝柔らかい）」ことを学習する。その結果、逆ヘッセ行列 $H$ は、硬いバネ方向の移動を抑制し、柔らかいポテンシャル方向の移動を促進するように探索方向 $\mathbf{d}_k$ を変形させる。これにより、バネ定数の大きさに関わらず、効率的に経路を最適化することが可能になる。Sheppardらのベンチマークでも、L-BFGSはバネ定数 $k$ の増加に対して収束性がほとんど悪化しないことが示されている。

3.2 鞍点近傍での対称性の回復N#

EB力の非保存性は、ポテンシャル力の接線成分を除去し、バネ力の接線成分を加えるという射影操作に起因する。

しかし、最適化が進行し、イメージ列がMEPに収束していくと、経路の接線ベクトル $\hat{\boldsymbol{\tau}}$ はポテンシャル曲面の等高線に対して直交するようになる（MEPの定義より）。この極限において、ポテンシャル力 $-\nabla V$ は元々経路に垂直な成分しか持たなくなり、射影操作による歪みが減少する。Sheppardらは、MEP近傍においてはNEB力のヤコビアン行列が対称行列に近づく（エルミート性が回復する）ため、BFGSによるヘッセ行列の近似が正当化されると考察している。

3.3 正定値性の強制による安定化#

BFGSの更新公式は、生成される行列 $H_k$ が常に正定値（Positive Definite）であることを保証する性質を持つ。正定値行列による変換は、勾配との内積が常に負（ $-\mathbf{F}^T H \mathbf{F} < 0$ ）となる「降下方向」を生成することを意味する。たとえ真のヤコビアンが非対称であったり、負の固有値（鞍点付近での曲率）を持っていたとしても、L-BFGSは正定値性を強制した「安定化されたモデル」を構築し続ける。これにより、MD法で見られるようなエネルギー保存に起因する振動や、探索方向の迷走を防ぎ、着実に鞍点へと収束させることができる。4. 他手法との比較とベンチマークSheppardらは、Global L-BFGSの性能を検証するため、以下の従来手法と比較を行った。

最急降下法 (SD): 実装は容易だが、硬いバネに対して極めて脆弱。

Quick-min (QM): MDベース。速度の射影により振動を抑制するが、曲率情報を利用しないため収束は一次。

FIRE: QMの改良版。構造最適化では強力だが、NEBのような拘束付き問題ではL-BFGSに劣る。

共役勾配法 (CG): 理論上は二次収束だが、非保存力場においては共役性の喪失により頻繁なリスタートが必要となり、効率が低下する。

ベンチマーク結果のハイライト#

LEPSポテンシャル（2次元モデル）は、他の全ての手法と比較して、力の評価回数が数分の一で済んだ。特に、高精度な収束（力が $10^{-3}$ 以下）を求めた場合、MDベースの手法との差は桁違いに広がった。

Pt(100)表面上の原子交換拡散:7原子が動く多次元系においても、L-BFGSはFIREやCGに比べて2倍以上の速さで収束した。初期経路が直線補間でエネルギーが高い状態からの緩和であっても、L-BFGSは安定して動作した。

Fe中の水素拡散:質量の異なる原子が混在する系（Stiffな振動モードとSoftな拡散モードの混在）においても、L-BFGSのPreconditioning効果が発揮され、高速な収束を示した。

結論#

Sheppard, Terrell, Henkelmanによる2008年の研究は、NEB法における最適化アルゴリズムの選択において、決定的な指針を与えた。結論として、Global L-BFGS法は、以下の理由によりNEB計算における最適な選択肢である。超一次収束性: ヘッセ行列の近似利用により、一次収束のMDベース手法（Quick-min, FIRE）を圧倒する収束速度を実現する。パラメータ非依存性: バネ定数による系の剛性（Stiffness）を自動的に補正するため、ユーザーによるパラメータ調整の負担が少ない。実装の普遍性: 全自由度を一括してベクトル化するアプローチにより、既存のL-BFGSライブラリをそのまま利用できる。この知見は、VASP、VTST Tools、ASE (Atomic Simulation Environment)、Quantum ESPRESSOなどの主要な第一原理計算コードにおけるNEB実装のデファクトスタンダードとなっており、今日における反応経路探索の高速化と高精度化を支える基盤技術となっている。

参考文献#

原著論文. Sheppard, R. Terrell, and G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys., 128, 134106 (2008).

Original NEB. Jónsson, G. Mills, and K. W. Jacobsen, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. J. Berne et al. (World Scientific, Singapore, 1998), p. 385.

L-BFGS Method. Nocedal, Updating Quasi-Newton Matrices with Limited Storage, Math. Comput., 35, 773 (1980).

FIRE Algorithm. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural Relaxation Made Simple, Phys. Rev. Lett., 97, 170201 (2006).