冗長内部座標系における自動化された鞍点探索アルゴリズムの数理と実装：Sellaについて

4234 words

21 minutes

冗長内部座標系における自動化された鞍点探索アルゴリズムの数理と実装：Sellaについて

2026-01-25

Computational Chemistry

Geometry Optimization

Saddle Point

Sella

Algorithm

Python

last_modified: 2026-01-25

生成AIによる自動生成記事に関する免責事項: 本記事は、Eric D. Hermes, Khachik Sargsyan, Habib N. Najm, Judit Zádorらによる原著論文 “Sella, an open-source automation-friendly molecular saddle point optimizer” (2022) の内容に基づき、大規模言語モデルによって作成された解説記事です。記述内容は提示された文献の範囲内において正確性を期していますが、厳密な証明や詳細な実装については、原著論文および関連する計算化学の専門書を参照してください。

1. 序論：反応経路探索における自動化の課題と座標系の選択#

化学反応における反応速度定数の算出や反応機構の解明において、ポテンシャルエネルギー曲面（PES）上の一次鞍点（first-order saddle point）、すなわち遷移状態（Transition State, TS）の特定は不可欠な手続きである。近年のエクサスケール・コンピューティングの到来に伴い、KinBotなどの自動化された反応経路探索フレームワークの重要性が増しているが、その中核となる「鞍点最適化」の堅牢性（robustness）がボトルネックとなっていた。

1.1 座標系の選択：デカルト座標 vs 内部座標#

構造最適化における座標系の選択は、収束効率に決定的な影響を与える。

デカルト座標: 実装が容易であり、原子間の接続関係（トポロジー）に依存しないため、凝縮相やクラスターのような密結合システムには適している。しかし、原子間の自由度が強く結合（カップリング）している分子系においては、PESの曲率が座標間で大きく異なるため、最適化効率が著しく低下する。
内部座標（Internal Coordinates）: 結合距離、結合角、二面角を用いる座標系。分子の化学的性質を反映しており、ポテンシャルエネルギーの二次微分（ヘシアン）の非対角項が小さくなる傾向があるため、分子系においてはデカルト座標よりも優れた収束性を示す。

しかし、内部座標を用いた自動化には重大な技術的障壁が存在した。それは、**「直線に近い結合角（near-linear angles）」の処理と、「トポロジーの変化」**への対応である。結合角が180度（または0度）に近づくと、その微分が定義不能あるいは数値的に不安定となり、計算の破綻を招く。従来、この問題は人間が手動でZ-matrixを修正するか、ダミー原子を追加することで対処されてきたが、これは完全自動化を目指すハイスループット計算においては許容されない介入である。

本稿で解説するアルゴリズム「Sella」は、冗長内部座標（Redundant Internal Coordinates）をベースとし、これらの特異点問題を数理的に回避・自動処理するメカニズムを実装したものである。

2. 数理的枠組み：座標変換と部分空間の射影#

Sellaのアルゴリズムは、複数の座標空間（デカルト空間、冗長内部空間、非冗長内部空間、制約付き空間）の間の射影と変換によって構成される。

2.1 座標空間の定義#

分子系を構成する原子数を $n$ とするとき、各空間は以下のように定義される。

デカルト座標空間 ( $x$ ): $3n$ 次元。最適化における独立変数は常にここにある。
冗長内部座標空間 ( $q$ ): $m$ 次元（通常 $m \ge 3n-6$ ）。結合、アングル、二面角の全セットを含む。
非冗長内部座標空間 ( $\breve{q}$ ): 冗長性を排除した空間。
制約付き空間と非制約空間 ( $\tilde{q}$ ): 制約条件を課した場合の部分空間。

2.2 射影行列の構築#

冗長内部座標から非冗長空間への変換には、ウィルソンB行列（ $J = dq/dx$ ）の特異値分解（SVD）を利用する。

J = \frac{dq}{dx} = \begin{bmatrix} N & R \end{bmatrix} \begin{bmatrix} T & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} Z^T \\ {Z'}^T \end{bmatrix} \quad (1)

ここで、 $T$ は非ゼロ特異値の対角行列、 $N$ は非冗長空間を張る正規直交行列、 $R$ は冗長空間（ヌル空間）を張る行列である。特異値が $10^{-6}$ 未満のものはゼロとみなされ、冗長空間へ割り当てられる。

この分解により、任意のベクトル $v$ および行列 $M$ は、冗長空間と非冗長空間の間で以下のように変換される。

\breve{v} = N^T v, \quad v = N \breve{v} \quad (2)

\breve{M} = N^T M N, \quad M = N \breve{M} N^T \quad (3)

2.3 制約条件の処理（Null-Space SQP）#

制約条件 $c(q) = 0$ が存在する場合、さらに空間を分割する必要がある。制約の数を $l$ とし、独立な制約の数を $p$ とするとき、制約空間への射影行列 $\tilde{P}$ と、最適化を行う非制約空間への射影行列 $\tilde{Q}$ は、以下のヤコビアンのSVDから得られる。

\frac{dc}{d\breve{q}} = \frac{dc}{dx} Z T^{-1} = \begin{bmatrix} \tilde{U} & \tilde{U}' \end{bmatrix} \begin{bmatrix} \tilde{W} & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} \tilde{P}^T \\ \tilde{Q}^T \end{bmatrix} \quad (5, 6)

これにより、全冗長空間 $q$ から、制約空間および非制約空間への直接射影行列 $P, Q$ が定義される。

P = N \tilde{P}, \quad Q = N \tilde{Q} \quad (7, 8)

この数学的定式化により、Sellaは制約条件（結合距離固定など）と最適化変数を厳密に直交分解し、数値的に安定した操作を可能にしている。

3. アルゴリズムのユニークな点：直線結合角とダミー原子の自動化#

本論文の貢献の一つに、最適化の途中であってもトポロジーの特異性を検知し、座標系を動的に再構築する機能である。

3.1 直線結合角の問題#

3原子A-B-Cがなす角 $\alpha_{abc}$ が $180^\circ$ または $0^\circ$ に近づくと、その微分係数は発散、あるいは不連続となる。また、この結合角を含む二面角は定義不能となる。従来手法では、この時点で計算が停止するか、デカルト座標へフォールバック（Q-Chem等の挙動）していた。

3.2 ダミー原子の自動挿入アルゴリズム#

Sellaは、結合角が $165^\circ$ を超えた（あるいは $15^\circ$ 未満になった）時点で「直線に近い」と判定し、以下のロジックで対処する。

不適切な二面角（Improper Dihedral）への置換: 中心原子BがA, C以外にも結合している原子Dを持つ場合、角度 $\alpha_{abc}$ を不適切な二面角 $\tau_{abdc}$ に置換する。
ダミー原子（Dummy Atom）の挿入: 原子BがA, Cとしか結合していない場合、ダミー原子 $X$ を系に追加する。

ダミー原子の配置ベクトル#

ダミー原子 $X$ の位置は、原子Bから距離 $1\mathrm{\AA}$ の位置に配置されるが、その方向ベクトル $\hat{\delta}_{bx}$ の決定には幾何学的な配慮がなされる。

通常ケース: 3原子が完全な直線でない場合、その平面に対する法線方向を用いる。 $\hat{\delta}_{bx} \parallel \hat{\delta}_{ba} \times \hat{\delta}_{bc}$
完全直線ケース（外積がゼロに近い場合）: 原子A-C間のベクトル $\hat{\delta}_{ac}$ と直交するデカルト軸 $\hat{w}$ を用いて射影する。 $\text{Direction} \propto (I - \hat{\delta}_{ac}\hat{\delta}_{ac}^T)\hat{w}$

自由度の相殺（Constraints）#

ダミー原子の導入は $3$ 自由度（デカルト座標成分）を増加させる。これによる物理的な意味の変化を防ぐため、Sellaは自動的に以下の $3$ つの制約を課す。

結合距離 $\delta_{bx} = 1.0 \mathrm{\AA}$
結合角 $\alpha_{abx} = 90^\circ$
結合角 $\alpha_{cbx} = 90^\circ$

これにより、本来の自由度は保存され、かつ特異点であった直線角 $\alpha_{abc}$ は、ダミー原子を介した不適切な二面角 $\tau_{abxc}$ として安定に記述される。

4. 反復的ヘシアン対角化と近似ヘシアンの構築#

鞍点探索においては、ヘシアン行列（エネルギーの二次微分）の最低固有値（負の曲率）方向へ進むことが求められる。しかし、 $3n \times 3n$ のヘシアンを完全に計算・対角化するコストは $O(N^4)$ 以上となることが多く、大規模系では現実的ではない。

Sellaは、完全なヘシアンを求めず、必要な固有ベクトルのみを抽出する**反復対角化（Iterative Diagonalization）**を採用している。

4.1 Olsen法による部分対角化#

ラグランジュ関数のヘシアン $H_{\mathcal{L}}$ の最小固有値と固有ベクトルを求めるために、一般化Davidson法の一種であるOlsen法を用いる。

具体的には、試行ベクトル空間 $\tilde{S}_i$ を拡張する際、以下の線形方程式を解いて修正ベクトル $\tilde{t}_i$ を求める。

(\tilde{I} - \tilde{z}_i \tilde{z}_i^T)(\tilde{B}_{\mathcal{L}} - \theta_i \tilde{I})(\tilde{I} - \tilde{z}_i \tilde{z}_i^T)\tilde{t}_i = -\tilde{r}_i \quad (14)

ここで $\tilde{B}_{\mathcal{L}}$ は現在の近似ヘシアン、 $\theta_i, \tilde{z}_i$ は現在の近似固有値・固有ベクトル、 $\tilde{r}_i$ は残差ベクトルである。この手法により、完全なヘシアンを構築することなく、ヘシアン・ベクトル積（Hessian-vector product）の評価のみで反応座標方向を特定できる。

4.2 数値微分の効率化#

ヘシアン・ベクトル積は、有限差分法を用いて計算される。

H s \approx \frac{\nabla \mathcal{E}(x + \eta s) - \nabla \mathcal{E}(x)}{\eta} \quad (12)

この操作は勾配計算1回分に相当し、完全なヘシアン（ $3n$ 回の勾配計算）に比べて圧倒的に低コストである。

4.3 マルチセカント法によるヘシアン更新#

対角化プロセスで得られた曲率情報（ $Y_{\mathcal{E}}$ と $S$ ）は、無駄なく近似ヘシアン $B$ の更新に利用される。ここでは、Multi-secant TS-BFGS更新式が採用されている。

B_{k+1} = B_k + U(Y_{\mathcal{E}} - B_k S)^T + (Y_{\mathcal{E}} - B_k S)U^T - U(Y_{\mathcal{E}} - B_k S)^T S U^T \quad (17)

これにより、探索が進むにつれて近似ヘシアンの精度が向上し、最適化ステップの質が高まる。

5. 最適化アルゴリズム：Null-Space SQPとRS-PRFO#

Sellaの最適化ステップは、制約条件を満たしつつ鞍点へ向かうため、Null-Space SQP（Sequential Quadratic Programming）とRS-PRFO（Restricted Step Partitioned Rational Function Optimization）を組み合わせた手法をとる。

5.1 変位ベクトルの分解#

全変位ベクトル $s$ は、制約多様体へ戻るためのステップ $s_P$ と、接空間内での最適化ステップ $s_Q$ に分解される。

s = s_P + s_Q \quad (23)

制約補正ステップ ( $s_P$ ): 制約 $c(q)=0$ を満たすための最小二乗解として求められる。
$\left(\frac{\partial c}{\partial q}\right) s_P = -c(q) \quad (24)$
最適化ステップ ( $s_Q$ ): 制約多様体の接空間（非制約空間）において、以下の外挿された勾配 $\tilde{g}_{ex}$ を用いて計算される。
$\tilde{g}_{ex} = \tilde{g} + Q^T B s_P \quad (26)$
この外挿により、制約を満たす過程で生じるポテンシャルエネルギーの変化を一次近似で予測し、振動を抑制する。

5.2 RS-PRFOによるステップ決定#

鞍点探索では、反応座標方向にはエネルギー最大化、それ以外の方向にはエネルギー最小化を行う必要がある。RS-PRFO法は、ヘシアンの固有モードに基づいて空間を「最大化部分空間」と「最小化部分空間」に分割し、それぞれに対して有理関数近似によるステップ制限（Trust Region）を設ける。

最大化ステップおよび最小化ステップは、それぞれ以下の固有値問題の解として導出される。

\begin{pmatrix} \alpha^2 \tilde{B}_{\mathcal{L}}^{(\text{sub})} & \alpha \tilde{g}_{ex}^{(\text{sub})} \\ \alpha \tilde{g}_{ex}^{(\text{sub})T} & 0 \end{pmatrix} \begin{pmatrix} \alpha^{-1} \tilde{s}_Q^{(\text{sub})} \\ 1 \end{pmatrix} = 2\mu \begin{pmatrix} \alpha^{-1} \tilde{s}_Q^{(\text{sub})} \\ 1 \end{pmatrix}

ここで $\alpha$ は信頼半径（Trust Region size）を制御するパラメータであり、無限大ノルム（Infinity-norm, $||s||_\infty$ ）に基づいて決定される。内部座標系においては、自由度の数が増減するため、ユークリッドノルム（2-norm）よりも無限大ノルムの方がスケーリング特性が良いとされる。

6. 実証結果とベンチマーク#

提案手法の有効性を検証するため、KinBotによって生成された500個の小規模有機分子の遷移状態構造候補を用いたベンチマークテストが行われた。比較対象として、NWChem、Q-Chem、Psi4の各主要量子化学パッケージの標準的な鞍点探索アルゴリズムが用いられた。

6.1 収束性能の比較#

図5および図6（原著論文参照）に示された結果は、以下の傾向を示している。

勾配評価回数の削減: Sella（デフォルト設定）は、完全なヘシアンを計算する他のコードと比較して、収束に至るまでの実効的な勾配評価回数（Effective Gradient Evaluations）が少ない。これは反復対角化によって、必要最小限の曲率情報のみを取得しているためである。
頑健性（Robustness）: 直線結合角を含む「病的な」構造（全体の約8.8%）に対しても、Sellaは座標系の自動再構築により計算を完遂した。特に、最適化の途中で結合角が直線化したケース（図中赤点）においても、Q-Chemのようなデカルト座標へのフォールバックなしに収束している。
コストの逆転: 非常に小規模な系や、初期構造が既に鞍点に近い場合、完全対角化を行う手法（またはSellaで厳密な収束基準を用いた場合）の方がステップ数は少なくなる傾向がある。しかし、系が大きくなるにつれ、ヘシアン計算コストが支配的になるため、部分対角化の優位性が増す。

6.2 制約付きスキャン（Dihedral Scan）#

SellaのNull-Space SQP実装は、反応経路の追跡（Relaxed Scan）においても高い効率を示した。5-exo-tet環化反応の遷移状態における二面角スキャン（図7）では、ステップ幅を細かく（ $1^\circ$ ）した場合、1ステップあたりの平均勾配評価回数は1.61回まで低下した。これは、前のステップのヘシアン情報を利用した外挿（Equation 26）が極めて有効に機能していることを示唆しており、理論下限値（1回）に近い高効率である。

7. 実装と利用方法#

SellaはPythonで記述されており、Atomic Simulation Environment (ASE) のOptimizerとして設計されている。これにより、ASEがサポートするあらゆる量子化学計算コード（NWChem, Gaussian, VASP, Quantum ESPRESSO等）とシームレスに連携可能である。

7.1 基本的な使用例#

以下に、Sellaを用いてNWChem経由で鞍点最適化を行う最小限のPythonスクリプトを示す（原著論文 Listing 1に基づく）。

#!/usr/bin/env python3
from ase.io import read
from ase.calculators.nwchem import NWChem
from sella import Sella

# 初期構造の読み込み
atoms = read('input_structure.xyz')

# 計算エンジンの設定（例：NWChem, HF/3-21G）
atoms.calc = NWChem()

# Sellaオプティマイザの初期化
# internal=True で内部座標の使用を明示
opt = Sella(atoms, internal=True)

# 最適化実行（収束条件: 最大力 < 0.01 eV/A）
opt.run(fmax=0.01)

このコードは、内部座標の生成、特異点のチェック、ヘシアンの反復対角化、および最適化ステップの実行をすべて自動的に処理する。

8. 結論#

本研究で提示されたSellaアルゴリズムは、内部座標を用いた鞍点探索における「自動化のラストワンマイル」であった特異点処理の問題を、ダミー原子の動的挿入とNull-Space SQPによる制約管理によって解決した。数学的な完全性: 冗長内部座標、非冗長空間、制約空間の間の射影を厳密に定式化し、数値的な安定性を確保した。計算効率: 反復的ヘシアン対角化とMulti-secant更新により、 $O(N^4)$ の解析的ヘシアン計算を回避しつつ、正確な反応座標方向への誘導を実現した。実用性: Python/ASEベースのオープンソース実装により、既存の多くの電子状態計算コードに対して、高度な鞍点探索機能をアドオンとして提供可能とした。この成果は、数千〜数万の反応を扱うハイスループットな反応探索プロジェクト（例えば燃焼化学や大気化学の自動モデリング）において、計算資源の大幅な節約と人的介入コストの削減に寄与するものである。

参考文献#

Hermes, E. D., Sargsyan, K., Najm, H. N., & Zádor, J. (2022). Sella, an open-source automation-friendly molecular saddle point optimizer. ChemRxiv. DOI: 10.26434/chemrxiv-2022-44r17-v2
Hermes, E. D., et al. (2019). Accelerated Saddle Point Refinement through Full Exploitation of Partial Hessian Diagonalization. J. Chem. Theory Comput., 15, 6536-6549.
Hermes, E. D., et al. (2021). Geometry optimization speedup through a geodesic approach to internal coordinates. J. Chem. Phys., 155, 094105.
Baker, J. (1997). Constrained optimization in delocalized internal coordinates. J. Comput. Chem., 18, 1079-1095
.Olsen, J., et al. (1990). Passing the one-billion limit in full configuration-interaction (FCI) calculations. Chem. Phys. Lett., 169, 463-472.