Home
4982 words
25 minutes
Cerjan-Miller法による遷移状態探索:Hessianを用いたポテンシャル曲面の局所二次近似と『上り坂』探索アルゴリズムの数理

最終更新:2026-01-03

注意: この記事は生成AIによって自動生成されたものです。内容は Charles J. Cerjan and William H. Miller, “On finding transition states”, J. Chem. Phys. 1981, 75, 2800 に基づいています。正確な学術的情報は原著論文を参照してください。

序論:ポテンシャルエネルギー曲面における鞍点探索の課題#

化学反応動力学や反応機構の解明において、Born-Oppenheimerポテンシャルエネルギー曲面(PES)上の特徴的な点、すなわち極小点(反応物および生成物)と、それらを繋ぐ一次の鞍点(遷移状態、Transition State: TS)を特定することは重要な課題である。 1980年代初頭、量子化学計算におけるエネルギー勾配(第一微分)の解析的導出法が確立されつつあったが、遷移状態の探索は、エネルギー最小化(構造最適化)と比較してはるかに困難な問題であった。エネルギー最小化であれば、勾配ベクトルに沿って「下り坂」を進むことで(最急降下法や共役勾配法などにより)比較的容易に極小点へ到達できる。しかし、遷移状態探索においては、反応座標に相当する特定のモードに対してはエネルギー極大(上り坂)であり、かつ他の全ての直交するモードに対してはエネルギー極小(下り坂)であるような点を見つけ出さなければならない。

1981年、Charles J. CerjanとWilliam H. Millerは、Journal of Chemical Physics 誌において、Hessian(エネルギーの二次微分行列)を用いた局所的な二次近似に基づく、ロバストな遷移状態探索アルゴリズムを提案した。本稿では、後に現代の幾何構造最適化手法(RFO法など)の基礎となったこの「Cerjan-Miller法」について、その数学的定式化、物理的意味、および歴史的意義を詳細に解説する。


1. 歴史的背景と従来法の限界#

CerjanとMillerの研究以前にも、遷移状態を探索するための手法はいくつか提案されていたが、それらは主に以下の二つのアプローチに分類され、それぞれに重大な欠点を抱えていた。

1.1 勾配ノルム最小化(Gradient Norm Minimization)#

これは、ポテンシャルエネルギー V(x)V(x) の勾配ベクトル g(x)=V(x)g(x) = \nabla V(x) のノルム g(x)|g(x)|、あるいはその二乗 g(x)2|g(x)|^2 を最小化する手法である。遷移状態では勾配がゼロになるため、ノルムの最小化によって探索可能であるという発想に基づく。 しかし、McIverやKomornickiらによって指摘されたように、この手法には以下の問題点があった。

  • 収束の不確実性: g(x)=0|g(x)|=0 は極小点、極大点、鞍点のすべてで成立するため、得られた点が目的の遷移状態である保証がない。
  • 計算コスト: ノルムの最小化には、元のポテンシャルの二次微分だけでなく三次微分の情報(あるいはその近似)が必要となり、計算コストが増大する。
  • 初期値依存性: 遷移状態に非常に近い初期構造を与えなければ、全く関係のない定留点へ収束するか、あるいは収束しないことが多い。

1.2 エネルギー最小化法の転用#

単純なエネルギー最小化手法を用いる場合、鞍点付近ではポテンシャル曲面が複雑な形状をしているため、通常の最適化アルゴリズムでは数値的な不安定性が生じやすい。また、特定の座標(反応座標)を固定して他を最適化する方法(Coordinate Driving)も存在するが、これは適切な反応座標をあらかじめ知っている必要があり、複雑な多原子分子系では「化学的直感」に頼らざるを得ないという限界があった。

CerjanとMillerは、これらの問題を克服するため、「極小点から出発して、自動的に反応経路に沿って『上り坂(uphill)』を歩き、遷移状態へ到達する」体系的なアルゴリズムを構築することを目的とした。


2. 理論的枠組みと数学的導出#

Cerjan-Miller法の核心は、現在の点における局所的な二次近似(Taylor展開)と、ステップサイズを制御するための制約条件付き最適化(Lagrange未定乗数法)の組み合わせにある。

2.1 ポテンシャルエネルギーの局所二次近似#

NN 個の原子からなる系(3N63N-6 自由度)において、現在の原子核座標を aa とする。ここから変位 Δx\Delta x だけ動いた点 x=a+Δxx = a + \Delta x におけるポテンシャルエネルギー V(x)V(x) を、点 aa 周りで二次までTaylor展開する。

V(a+Δx)V(a)+DTΔx+12ΔxTKΔxV(a + \Delta x) \approx V(a) + \mathbf{D}^T \Delta x + \frac{1}{2} \Delta x^T \mathbf{K} \Delta x

ここで、各項は以下のように定義される。

  • V(a)=V0V(a) = V_0: 現在の点におけるエネルギー
  • D=V(a)\mathbf{D} = \nabla V(a): エネルギー勾配ベクトル(Gradient, 第一微分)
  • K=2V(a)\mathbf{K} = \nabla^2 V(a): 力定数行列(Hessian, 第二微分)

2.2 ステップ方向の定性的考察#

もしポテンシャル曲面を一次近似(線形近似)のみで考えると、V(x)V0+DTΔxV(x) \approx V_0 + \mathbf{D}^T \Delta x となり、等高線は直線群となる。この場合、最も急な上り坂は勾配ベクトル D\mathbf{D} の方向である。しかし、遷移状態探索において単純に勾配ベクトルに沿って登ることは、一般に失敗することが知られている。これは、反応経路(谷底の最深部を繋ぐ線)が曲がっている場合、勾配方向へ進むと谷の壁を登ってしまい、反応経路から外れてしまうためである(いわゆる「defocusing」現象)。

Hessian K\mathbf{K} を含めた二次近似を導入することで、ポテンシャル曲面の曲率(谷の形状)を考慮することが可能になる。CerjanとMillerは、固定されたステップ長 Δx=Δ|\Delta x| = \Delta という制約の下で、局所的な二次モデル上の定留点を探索する方針を採った。

2.3 Lagrange未定乗数法による定式化#

ステップ長を Δ\Delta に固定するという制約条件 ΔxTΔx=Δ2\Delta x^T \Delta x = \Delta^2 を課した上で、二次近似されたエネルギーの極値を求める。これは以下のLagrange関数 L(Δx,λ)L(\Delta x, \lambda) の停留点を求める問題に帰着する。

L(Δx,λ)=V0+DTΔx+12ΔxTKΔxλ2(ΔxTΔxΔ2)L(\Delta x, \lambda) = V_0 + \mathbf{D}^T \Delta x + \frac{1}{2} \Delta x^T \mathbf{K} \Delta x - \frac{\lambda}{2} (\Delta x^T \Delta x - \Delta^2)

ここで λ\lambda はLagrange未定乗数である。停留条件は、Δx\Delta x および λ\lambda による偏微分がゼロになることである。

LΔx=D+KΔxλΔx=0(1)\frac{\partial L}{\partial \Delta x} = \mathbf{D} + \mathbf{K} \Delta x - \lambda \Delta x = 0 \quad \cdots (1) Lλ=12(Δ2ΔxTΔx)=0(2)\frac{\partial L}{\partial \lambda} = \frac{1}{2} (\Delta^2 - \Delta x^T \Delta x) = 0 \quad \cdots (2)

式(1)を変形すると、以下の固有値方程式のような形式が得られる。

(KλI)Δx=D(\mathbf{K} - \lambda \mathbf{I}) \Delta x = - \mathbf{D}

したがって、変位ベクトル Δx\Delta x は形式的に以下のように表される。

Δx(λ)=(KλI)1D=(λIK)1D(3)\Delta x(\lambda) = - (\mathbf{K} - \lambda \mathbf{I})^{-1} \mathbf{D} = (\lambda \mathbf{I} - \mathbf{K})^{-1} \mathbf{D} \quad \cdots (3)

この式は、現代の「信頼半径法(Trust Region Method)」や「シフト付きNewton法(Shifted Newton Method)」の基本式そのものである。λ=0\lambda = 0 の場合、これは通常のNewton-Raphsonステップ Δx=K1D\Delta x = - \mathbf{K}^{-1} \mathbf{D} に一致する。

2.4 スペクトル分解による解析#

Hessian行列 K\mathbf{K} は実対称行列であるため、直交行列 U\mathbf{U} によって対角化可能である。

UTKU=k=diag(k1,k2,,kF)\mathbf{U}^T \mathbf{K} \mathbf{U} = \mathbf{k} = \text{diag}(k_1, k_2, \dots, k_F)

ここで kik_i はHessianの固有値、FF は自由度である。また、勾配ベクトル D\mathbf{D} をHessianの固有ベクトル空間へ射影したものを d=UTD\mathbf{d} = \mathbf{U}^T \mathbf{D} とすると、式(3)のステップベクトルの二乗ノルム Δ2\Delta^2 は、λ\lambda の関数として以下の有理関数(Rational Function)で表される。

Δ2(λ)=Δx2=DT(λIK)2D=i=1Fdi2(λki)2(4)\Delta^2(\lambda) = |\Delta x|^2 = \mathbf{D}^T (\lambda \mathbf{I} - \mathbf{K})^{-2} \mathbf{D} = \sum_{i=1}^{F} \frac{d_i^2}{(\lambda - k_i)^2} \quad \cdots (4)

同様に、このステップを採用した場合のエネルギー変化 ΔV(λ)=V(Δx)V0\Delta V(\lambda) = V(\Delta x) - V_0 は以下のようになる。

ΔV(λ)=i=1Fdi2(λ12ki)(λki)2(5)\Delta V(\lambda) = \sum_{i=1}^{F} \frac{d_i^2 (\lambda - \frac{1}{2}k_i)}{(\lambda - k_i)^2} \quad \cdots (5)

2.5 未定乗数 λ\lambda の選択と「上り坂」探索#

式(4)および式(5)の挙動を解析することが、Cerjan-Miller法の核心部分である。 関数 Δ2(λ)\Delta^2(\lambda) は、λ\lambda がHessianの固有値 kik_i に近づくと発散する特異点を持つ。したがって、Δ2(λ)\Delta^2(\lambda) のグラフは、固有値 kik_i によって区切られた複数の分枝(branch)を持つことになる。

CerjanとMillerは、極小点から遷移状態へ向かうための適切なステップを得るために、λ\lambda をどのように選ぶべきかを以下のように論じた。

A. 極小点近傍(全ての ki>0k_i > 0 の場合)#

極小点付近ではHessianの全ての固有値は正である。このとき、関数 Δ2(λ)\Delta^2(\lambda) は区間 (ki,ki+1)(k_i, k_{i+1}) の間で極小値を持つ。 論文中のFig. 3に示されているように、Δ2(λ)\Delta^2(\lambda) の極小点に対応する λ\lambda の値(これを λ0\lambda_0 とする)を選択することが推奨される。 特に、最も低い固有値 k1k_1 と次の固有値 k2k_2 の間にある λ0\lambda_0 を選ぶと、以下の性質が得られることが示される。

  1. Hessianの最低モード方向への寄与: 式(3)の分母 (λki)(\lambda - k_i) を見ると、λ\lambdak1k_1 に近い正の値であれば、最低固有値 k1k_1 に対応するモード(最も曲率が緩やかな方向=反応経路の谷底方向)の成分が支配的となる。
  2. 上り坂への誘導: 式(5)のエネルギー変化において、適切な λ0\lambda_0 を選ぶことで ΔV>0\Delta V > 0、すなわちエネルギーが増加する方向へのステップが生成される。これは、極小点から脱出しようとする動きに対応する。

B. 遷移状態近傍(最低固有値 k1<0k_1 < 0 の場合)#

探索が進み、遷移状態に近づくと、反応座標に対応する最低固有値 k1k_1 が負になる(虚振動数を持つ)。 この領域(鞍点領域)に入ると、もはや「上り坂」を歩く必要はなく、速やかに鞍点へ収束させることが望ましい。 k1<0k_1 < 0 の場合、Newton-Raphson法に対応する λ=0\lambda = 0 が適切な選択肢となる。式(3)において λ=0\lambda=0 とすると、負の固有値 k1k_1 に対応する成分は、

Δx1=d1k1=d1k1\Delta x_1 = \frac{d_1}{-k_1} = \frac{d_1}{|k_1|}

となり、これは勾配 d1d_1 の方向(エネルギーが高い方向)へ進むことを意味する。一方で、正の固有値 ki(i>1)k_i (i > 1) に対応する成分はエネルギーを最小化する方向へ働く。これにより、反応座標については極大化、その他の座標については極小化が行われ、自然に鞍点へと収束する。

2.6 アルゴリズムの要約#

Cerjan-Millerのアルゴリズムは以下のように要約できる。

  1. 現在の点 aa での勾配 D\mathbf{D} とHessian K\mathbf{K} を計算する。
  2. Hessianを対角化し、固有値 kik_i と固有ベクトルを求める。
  3. 方程式 ddλΔ2(λ)=2di2(λki)3=0\frac{d}{d\lambda} \Delta^2(\lambda) = -2 \sum \frac{d_i^2}{(\lambda - k_i)^3} = 0 の解(根)を探す。
    • 初期段階(極小点近傍)では、λ>k1\lambda > k_1 となる最小の根 λ0\lambda_0 を採用する。
    • もし λ0\lambda_0 が負になる、あるいは k1k_1 が負になった場合は、状況に応じて λ=0\lambda=0 (Newton-Raphsonステップ)あるいは適切な負の固有値に対応する分枝へ切り替える。
  4. 決定された λ\lambda を用いてステップ Δx\Delta x を計算し、座標を更新する。
  5. 収束基準(勾配ノルムが十分小さい等)を満たすまで反復する。

3. 実利的な適用例と成果#

論文では、提案手法の有効性を検証するために、モデルポテンシャルおよび実際の分子系への適用結果が示されている。

3.1 2次元モデルポテンシャル(Muller-Brown型)#

数式 V(x,y)=(aby2)x2ex2+c2y2V(x,y) = (a - by^2)x^2 e^{-x^2} + \frac{c}{2}y^2 で表される2次元ポテンシャルを用いたテストが行われた。このポテンシャルは、原点 (0,0)(0,0) に極小値を持ち、(±1,0)(\pm 1, 0) 付近に遷移状態を持つ。 特筆すべきは、初期ステップを意図的に「間違った方向(y軸方向)」へ踏み出した場合でも、アルゴリズムが自動的にHessianの固有ベクトルに沿って方向修正を行い、正しくx軸方向にある遷移状態へ到達した点である(論文 Fig. 6)。これは、勾配のみに頼る手法では達成困難な挙動であり、Hessian情報の有効性を如実に示している。

3.2 ビニリデン-アセチレン転位反応#

H2C=C:HCCH\text{H}_2\text{C=C:} \leftrightarrow \text{HC}\equiv\text{CH} この1,2-水素移動反応は、有機化学的に興味深い多原子分子系への適用例である。 Cerjanらは、STO-2G基底系を用いた計算において、ビニリデン(極小点)から出発し、わずか8ステップで遷移状態構造へ到達することに成功した(Table I)。比較対象として挙げられているPoppingerの勾配ノルム最小化法では、遷移状態に近い適切な初期推定値を与えたとしても21回以上の勾配計算が必要であったとされ、Cerjan-Miller法の高い効率性が実証された。

3.3 ホルムアルデヒドの解離と異性化#

H2COH2+CO(解離)\text{H}_2\text{CO} \to \text{H}_2 + \text{CO} \quad (\text{解離}) H2COHCOH(異性化)\text{H}_2\text{CO} \to \text{HCOH} \quad (\text{異性化}) ホルムアルデヒドの系では、初期の摂動(perturbing step)の与え方や、追跡するHessianの固有モード(固有ベクトル)を選択することで、解離反応と異性化反応という異なる二つの遷移状態を同一の極小点から分別して探索できることが示された。 解離反応の遷移状態へは11ステップ、異性化反応へは9ステップで到達しており、Hessianの固有モードが反応経路の「谷」の方向を極めて正確に示唆していることがわかる。特に、初期段階では水素原子の角度変化(面外変角など)が大きく寄与し、その後結合距離が緩和するという物理的に妥当な反応経路が自動的に生成されている(Fig. 7, 8)。

3.4 反応経路(Intrinsic Reaction Coordinate)の追跡#

論文の後半(Fig. 9)では、遷移状態から極小点へ向かう「下り坂」の経路探索(IRC探索に相当)についても言及している。 最急降下法(単純な勾配法)では、谷底に対して振動(oscillation)しながら進むという数値的な不安定性がしばしば問題になる。Cerjan-Miller法の定式化(Hessianを用いたステップ)は、この振動を抑制し、滑らかに反応経路を追跡する上でも有効であることが、HNC→HCN異性化反応の例を用いて示された。これは、Hessianが谷の断面形状(放物線)を近似しているため、谷底の中心へ向かう補正が自然に組み込まれているからである。


4. 議論と歴史的意義#

4.1 「Gradient Extremal」との関係#

本論文の手法は、しばしば「Gradient Extremal Following (GEF)」の文脈で語られるが、厳密にはGEFの定義(勾配ノルムが等高線上で極値をとる軌跡)そのものを追跡するアルゴリズムではない。 しかし、Hessianの最小固有値に対応する固有ベクトルに沿って進むという戦略は、結果として反応経路(Minimum Energy Path, MEP)やGradient Extremalに近い軌跡を描くことが多い。Ruedenbergらが提唱した厳密なGEF理論と、Cerjan-Millerの実用的なアルゴリズムは、80年代から90年代にかけての反応経路探索手法の発展において、車の両輪のような役割を果たした。

4.2 現代的手法への系譜:RFO法へ#

CerjanとMillerが導入した Δx=(λIK)1D\Delta x = (\lambda \mathbf{I} - \mathbf{K})^{-1} \mathbf{D} という式と、λ\lambda を有理関数 Δ2(λ)\Delta^2(\lambda) の解析から決定する手法は、後にSimonsらによる Rational Function Optimization (RFO) 法へと直接的に発展した。 RFO法は、Lagrange関数をより一般的なPadé近似の形式として解釈し直したものであり、現在の主要な量子化学計算パッケージ(Gaussian, GAMESS, ORCAなど)における構造最適化および遷移状態探索の標準的なアルゴリズムとして採用されている。

また、ステップサイズを制御する Δ2\Delta^2 の制約は、現代の最適化理論における 信頼領域(Trust Region) の概念と等価である。Hessianモデルが信頼できる範囲内(球状領域)でステップを制限することで、大域的な収束安定性を劇的に向上させた点は、計算化学のみならず数理最適化の分野から見ても先駆的な成果であった。

4.3 計算コストの変遷#

論文中では、Hessian行列の計算コストが高いこと(当時は有限差分で求めていたため、多数の勾配計算が必要であった)が懸念点として挙げられていた。しかし、著者らは「近い将来、解析的二次微分(Analytic Hessian)の計算手法が普及すれば、この手法が真価を発揮する」と予見していた。 この予見は完全に的中した。Popleらによる解析的二次微分法の確立とコンピュータの性能向上により、Hessianを各ステップで計算(あるいは更新)することは日常的な処理となり、Cerjan-Miller法のアプローチは現代計算化学のデファクトスタンダードとなったのである。


結論#

CerjanとMillerによる1981年の論文は、単なる一つの探索アルゴリズムの提案にとどまらず、ポテンシャルエネルギー曲面上の探索という問題に対して、「Hessianを用いた局所情報の最大活用」と「信頼半径によるステップ制御」という、現代に通じる二つの基本指針を確立した記念碑的研究である。 彼らが導出したLagrange未定乗数法に基づくステップ決定式は、極小点からの「上り坂」探索を可能にし、化学者の直感に依存しない自動的な遷移状態探索への道を開いた。このアルゴリズムの数理的基盤は、半世紀近く経過した現在もなお、最新の計算化学ソフトウェアの深層で稼働し続けており、その実利的な価値は計り知れない。


参考文献#

  1. C. J. Cerjan and W. H. Miller, “On finding transition states”, Journal of Chemical Physics, 75, 2800–2806 (1981).
  2. J. W. McIver, Jr. and A. Komornicki, “Structure of transition states in organic reactions. General theory and an application to the cyclobutene-butadiene isomerization”, Journal of the American Chemical Society, 94, 2625–2633 (1972).
  3. D. Poppinger, “Geometry optimization in ab initio molecular orbital theory”, Chemical Physics Letters, 35, 550–554 (1975).
  4. J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Binkley, “Derivative studies in hartree-fock and moller-plesset theories”, International Journal of Quantum Chemistry Symposium, 13, 225–241 (1979).
  5. K. Müller and L. D. Brown, “Location of Saddle Points. I. Method and Application to H2 + H —> H + H2”, Theoretica Chimica Acta, 53, 75–93 (1979).
  6. W. H. Miller, N. C. Handy, and J. E. Adams, “Reaction path Hamiltonian for polyatomic molecules”, Journal of Chemical Physics, 72, 99–112 (1980).
Cerjan-Miller法による遷移状態探索:Hessianを用いたポテンシャル曲面の局所二次近似と『上り坂』探索アルゴリズムの数理
https://ss0832.github.io/posts/20260103_compchem_gef/
Author
ss0832
Published at
2026-01-03

Related Posts

ポテンシャル曲面交差上の最小エネルギー点を特定するための直接法:射影演算子を用いた制約なし最適化アルゴリズム
2026-01-04
M. J. Bearpark, M. A. Robb, H. B. Schlegelによる論文 *Chem. Phys. Lett.* **1994**, *223*, 269-274 の包括的解説。光化学反応の機構解明に不可欠な円錐交差(Conical Intersection)および項間交差の最低エネルギー点を、Lagrange未定乗数法を用いずに探索する「直接法(Direct Method)」の理論的導出、アルゴリズムの実装、およびベンチマーク結果について詳述する。
計算化学における有効へシアンと部分へシアンの数理:シューア補行列による統一的理解と使い分け
2026-01-06
巨大分子系の一部(部分系)を解析する際、周囲の環境を「緩和させる(有効へシアン)」か「固定する(部分へシアン)」かは、解析目的によって使い分けることが望ましい。本稿では、有効へシアンが数学的に「シューア補行列」と対応していることを示した上で、反応性解析には有効へシアンが、局所電子応答の評価には部分へシアンが適していると考えられる物理的・数理的根拠を整理する。
同一対称性を持つ2つのポテンシャルエネルギー曲面の交差:ラグランジュ未定乗数法を用いた体系的特性化と最小エネルギー交差探索
2026-01-04
M. Riad ManaaとDavid R. Yarkonyによる論文 *J. Chem. Phys.* **1993**, *99*, 5251 の詳細な解説。同一の空間・スピン対称性を持つ2つの電子状態間の円錐交差(Conical Intersection)を、ラグランジュ未定乗数法を用いて直接探索するアルゴリズムの数学的定式化、Hessian行列を用いたニュートン・ラフソン法の実装、およびHCOラジカルを用いたベンチマーク結果について、学術的な観点から深掘りする。
分岐平面更新法(UBP法)による最小エネルギー円錐交差の探索:非断熱結合ベクトルを用いない効率的アルゴリズムの理論と実装
2026-01-04
Satoshi Maeda, Koichi Ohno, Keiji Morokumaによる論文 *J. Chem. Theory Comput.* **2010**, *6*, 1538-1545 の詳細な学術的解説。最小エネルギー円錐交差(MECI)探索において計算コストのボトルネックとなる非断熱結合ベクトル(CDV)の算出を回避し、分岐平面の更新によって効率的な最適化を実現するUBP法の数学的定式化、歴史的背景、およびベンチマーク結果について詳述する。
Coordinate Driving法とVADERプログラムによる反応経路ネットワークの自動探索:理論的枠組みと実装
2026-01-03
Martin Černohorský, Jaroslav Kočaらによる1999年の論文 'VADER: New Software for Exploring Interconversions on Potential Energy Surfaces' に基づき、Single Coordinate Driving (SCD) 法の拡張と、それを実装したプログラムVADERについて詳説する。ポテンシャルエネルギー曲面(PES)のトポロジー解析、反応経路の自動探索アルゴリズム、およびホルムアルデヒドやtert-ブチルハライド系への適用事例を数理的・化学的観点から包括的に解説する。
ポテンシャルエネルギー局面上の2次最急降下法:局所二次近似 (LQA) におけるSun-Ruedenberg法の定式化と定量的評価
2026-01-03
1993年、Jun-Qiang SunとKlaus Ruedenbergによって提案された、ポテンシャルエネルギー局面(PES)上の最急降下経路(IRC)を決定するための新規2次アルゴリズムの包括的解説。本稿では、局所的な二次Taylor展開に基づく厳密な最急降下線の接続手法、および展開中心点(Expansion Center)、接続点(Junction Point)、予測点(Predicted Point)の3点を区別することによる精度向上の数学的基盤について詳述する。