解析的ヘシアンなしでの反応経路追跡：Sun-Ruedenbergによる二次最急降下法 (QSD) の実用アルゴリズム

2508 words

13 minutes

解析的ヘシアンなしでの反応経路追跡：Sun-Ruedenbergによる二次最急降下法 (QSD) の実用アルゴリズム

2026-01-03

Physical Chemistry

Computational Chemistry

last_modified: 2026-01-03

免責事項: 本記事は、Jun-Qiang Sun, Klaus Ruedenbergによる論文 The Journal of Chemical Physics 1993, 99, 5269-5275 を基に、生成AIによって作成された解説記事です。数式や論理の正確な理解のためには、必ず原著論文を参照してください。

1. 序論：LQA法の「ヘシアンの壁」とQSD法の提案#

1988年にPageとMcIverによって提案された局所二次近似 (LQA) 法は、反応経路（IRC）探索におけるステップ幅の制約を大幅に緩和し、硬い（stiffな）微分方程式に対するロバスト性を示した。しかし、LQA法の厳密な適用には、各ステップにおける解析的ヘシアン（二階微分行列） の計算が必要不可欠であった。

量子化学計算において、エネルギーの二階微分を解析的に求める計算コストは、エネルギーそのものや勾配（一階微分）の計算に比べて極めて高い（通常、原子数の4乗〜5乗のオーダーで増大する）。そのため、LQA法を巨大分子系に適用することは実質的に困難であった。

1993年、SunとRuedenbergは、一連の論文（Part I〜III）において、LQA法を数学的により厳密化した 二次最急降下 (Quadratic Steepest Descent; QSD) 法を提唱した。特に、本稿で解説する Part II の論文では、「解析的ヘシアンを用いずに、いかにしてQSD法を実行するか」 という実用上の最重要課題に対する解法が提示されている。

2. 理論的背景：QSD法の基礎方程式#

まず、Part Iで確立されたQSD法の基本原理を概観する。ポテンシャルエネルギー曲面 $E(\mathbf{x})$ を、展開中心 $\mathbf{x}_k$ の周りで二次形式近似する。

E(\mathbf{x}) \approx E_k + \mathbf{g}_k^T (\mathbf{x} - \mathbf{x}_k) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_k)^T \mathbf{H}_k (\mathbf{x} - \mathbf{x}_k) \quad (1)

ここで $\mathbf{g}_k$ は勾配ベクトル、 $\mathbf{H}_k$ はヘシアン行列である。 IRCは最急降下経路（ $\frac{d\mathbf{x}}{ds} = -\frac{\mathbf{g}}{|\mathbf{g}|}$ ）として定義される。Page-McIverのLQA法と同様に、QSD法でもこの近似曲面上での経路を解析的に求めるが、SunとRuedenbergは、次の展開点 $\mathbf{x}_{k+1}$ を決定するためのより厳密な幾何学的基準を導入した。

QSDステップにおける変位 $\Delta \mathbf{x}(s)$ は、ヘシアンの固有ベクトル $\mathbf{v}_i$ と固有値 $h_i$ を用いて以下のように表される。

\Delta \mathbf{x}(s) = \sum_{i} c_i(s) \mathbf{v}_i \quad (2)

係数 $c_i(s)$ は、Page-McIver法と同様に指数関数的な減衰（または成長）項を含むが、QSD法ではステップ間の接続における接線の連続性などがより深く考慮されている。

3. 具体的な求め方：更新法を用いたQSDアルゴリズム#

Part IIの核心は、毎ステップ $\mathbf{H}_k$ を真面目に計算する代わりに、ヘシアン更新法 (Hessian Update Method) と 数値微分 を組み合わせるハイブリッド戦略にある。

3.1 戦略の概要#

著者は、経路上の点を以下の2種類に分類した。

再生点 (Regeneration Points, $P_m$ ): ヘシアンを「完全に」計算し直す点。ただし、解析的微分ではなく、勾配を用いた数値微分によって精度良く求める。
中間点 (Intermediate Points, $Q_k$ ): 再生点の間にある点。ここではヘシアンを計算せず、前ステップのヘシアンと最新の勾配情報を用いて、行列を更新 (Update) する。

この戦略により、高い精度を維持しつつ、高コストなヘシアン計算の頻度を劇的に減らすことが可能となる。

3.2 ヘシアンの数値的再生 (Numerical Regeneration)#

再生点 $P_m$ において、解析的ヘシアンが利用できない場合、著者は勾配の差分による数値微分を採用した。 $N$ 個の原子を持つ系（自由度 $n=3N$ ）において、各座標軸方向に微小変位 $\delta$ を加えた点での勾配 $\mathbf{g}(\mathbf{x} + \delta \mathbf{e}_j)$ を計算し、次式でヘシアンの列ベクトルを求める。

(\mathbf{H})_{ij} \approx \frac{g_i(\mathbf{x} + \delta \mathbf{e}_j) - g_i(\mathbf{x} - \delta \mathbf{e}_j)}{2\delta} \quad (3)

この操作には $2n$ 回の勾配計算が必要となるが、毎ステップ行うわけではないため、トータルのコストは抑制される。

3.3 ヘシアン更新アルゴリズム (Update Schemes)#

中間点では、準ニュートン法などで用いられる更新公式を適用する。反応経路（特に遷移状態付近）では、ヘシアンが負の固有値を持つ（不定符号である）ため、通常の最小化で用いられるBFGS法などが常に最適とは限らない。論文では、以下の2つの更新式が検討されている。

MS (Murtagh-Sargent) / PSB (Powell-Symmetric-Broyden) 更新: 対称性を保ちつつ、不定符号のヘシアンにも対応可能。 $\mathbf{H}_{k+1} = \mathbf{H}_k + \frac{(\mathbf{y}_k - \mathbf{H}_k \mathbf{s}_k)(\mathbf{y}_k - \mathbf{H}_k \mathbf{s}_k)^T}{(\mathbf{y}_k - \mathbf{H}_k \mathbf{s}_k)^T \mathbf{s}_k} \quad (MS)$
BFGS / DFP 更新: 正定値性を保持する傾向があるため、極小点探索では強力だが、IRC探索では注意が必要。しかし、著者の検証では、ステップサイズが十分に小さければBFGSも有効に機能することが示された。

ここで、 $\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k$ は変位ベクトル、 $\mathbf{y}_k = \mathbf{g}_{k+1} - \mathbf{g}_k$ は勾配の差ベクトルである。

3.4 再生基準 (Regeneration Criterion)#

いつヘシアンを再生するかという「トリガー」の設計が、効率と精度のバランスを決定する。SunとRuedenbergは、直前の再生点からの累積経路長 $S$ に基づく基準を提案した。

S_{cum} = \sum |\Delta \mathbf{x}| > S_{max} \quad (4)

$S_{max}$ は「信頼半径」のような役割を果たし、更新されたヘシアンの精度が劣化する前に、正しいヘシアンを再計算するよう制御する。論文では、典型的な化学反応において $S_{max} \approx 0.1 \sim 0.5 \, \text{bohr}$ 程度が推奨されている。

4. アルゴリズムの実装フロー#

本論文で提案されたQSD法（Part II）の実装ロジックを、擬似コード形式で記述する。

def QSD_without_Analytic_Hessian(x_start, direction, S_max, tolerance):
    """
    Sun-Ruedenberg (1993) Part IIに基づくQSD法のアルゴリズム概要
    
    Args:
        x_start: 遷移状態（または開始点）の座標
        direction: 反応進行方向（遷移ベクトル）
        S_max: ヘシアン再生を行う累積距離の閾値
        tolerance: 収束判定基準
    """
    
    # 1. 初期化
    x_current = x_start
    g_current = calculate_gradient(x_current)
    
    # 初期ヘシアンは数値微分で正確に求める
    H_current = calculate_numerical_hessian(x_current)
    
    path_points = [x_current]
    accumulated_distance = 0.0
    
    while True:
        # 2. 収束判定 (勾配ノルムやエネルギー変化など)
        if check_convergence(g_current, tolerance):
            break
            
        # 3. QSDステップの計算 (Part Iの理論に基づく)
        # 局所二次近似曲面上での最急降下経路を解析的に解く
        # delta_x = Sum [ c_i(s) * v_i ]
        delta_x = compute_QSD_step(g_current, H_current)
        
        step_size = norm(delta_x)
        x_next = x_current + delta_x
        g_next = calculate_gradient(x_next)
        
        # 4. 累積距離の更新
        accumulated_distance += step_size
        
        # 5. ヘシアンの取り扱い (Part IIの核心)
        if accumulated_distance > S_max:
            # 基準距離を超えた場合：ヘシアンを数値的に「再生」する
            H_next = calculate_numerical_hessian(x_next)
            accumulated_distance = 0.0  # カウンタリセット
            method = "Regeneration"
        else:
            # 基準距離以内の場合：ヘシアンを「更新」する
            # s = x_next - x_current
            # y = g_next - g_current
            # H_next = Update(H_current, s, y) using MS or BFGS formula
            H_next = update_hessian(H_current, x_next - x_current, g_next - g_current)
            method = "Update"
            
        # 6. 次のステップへの準備
        x_current = x_next
        g_current = g_next
        H_current = H_next
        
        path_points.append(x_current)
        print(f"Step: {len(path_points)}, Dist: {step_size:.4f}, Method: {method}")

    return path_points

5. 実利的な成果と評価#

5.1 精度とコストのトレードオフ#

論文では、Müller-BrownポテンシャルやHCN $\to$ HNC異性化反応などのモデル系を用いて検証が行われた。

解析的ヘシアンQSD: 最も正確だが、コスト大。

更新QSD (本手法): 解析的QSDとほぼ同等の軌跡を描きつつ、勾配計算回数を大幅に削減することに成功した。

Gonzalez-Schlegel (GS) 法との比較: GS法（拘束付き最適化に基づく手法）と比較しても、QSD法は特に曲率が大きい領域（反応経路が急カーブする領域）において、より忠実に最急降下経路を追跡できることが示された。

5.2 現代計算化学への影響#

この「ヘシアン更新を用いたQSD法」は、その後の量子化学パッケージ（GAMESS, MOLPRO等）におけるIRC計算の標準的なオプションとして採用された。特に、以下の点が実務上重要である。遷移状態からの脱出: 遷移状態近傍では正確な曲率情報が必須であるため、最初は数値ヘシアンを用い、平坦な領域に入ったら更新法に切り替えるという柔軟な運用が可能になった。巨大系への適用: DFT計算など、解析的ヘシアンが利用可能であっても高コストな場合において、この更新アルゴリズムは計算時間を数分の一に短縮する効果がある。

6. 結論#

SunとRuedenbergによる1993年の研究（Part II）は、理想的な理論（QSD）を現実的な計算リソースで実行可能にするための「アルゴリズムの架け橋」を構築したものである。彼らは、「必要なときだけ高コストな計算（再生）を行い、それ以外は数学的な近似（更新）で済ませる」という戦略を、厳密な基準 $S_{max}$ に基づいて定式化した。これにより、IRC計算は単なる理論的概念の確認作業から、複雑な多原子分子の反応機構を日常的に解析するための実用的なツールへと進化したのである。

参考文献#

Sun, J. Q.; Ruedenberg, K. “Quadratic steepest descent on potential energy surfaces. I. Basic formalism and quantitative assessment”, J. Chem. Phys. 1993, 99, 5257-5268.
Sun, J. Q.; Ruedenberg, K. “Quadratic steepest descent on potential energy surfaces. II. Reaction path following without analytic Hessians”, J. Chem. Phys. 1993, 99, 5269-5275.
Page, M.; McIver, J. W., Jr. J. Chem. Phys. 1988, 88, 922.
Gonzalez, C.; Schlegel, H. B. J. Chem. Phys. 1989, 90, 2154.
Murtagh, B. A.; Sargent, R. W. H. Comput. J. 1970, 13, 185.