多参照配置間相互作用法とReduced Gradient Following法を組み合わせたエネルギー極小点および鞍点の決定：H2COのS0曲面とアセチレンのT1, T2曲面への適用

4840 words

24 minutes

多参照配置間相互作用法とReduced Gradient Following法を組み合わせたエネルギー極小点および鞍点の決定：H2COのS0曲面とアセチレンのT1, T2曲面への適用

2026-01-03

Theoretical Chemistry

Computational Chemistry

最終更新：2026-01-03

注意: この記事は生成AIによって自動生成されたものです。内容は Michal Dallos, Hans Lischka, Elizete Ventura do Monte, Michael Hirsch, Wolfgang Quapp, “Determination of Energy Minima and Saddle Points Using Multireference Configuration Interaction Methods in Combination with Reduced Gradient Following”, J. Comput. Chem. 2002, 23, 576-583 に基づいています。正確な学術的情報は原著論文を参照してください。

序論：多参照理論における停留点探索の課題と展望#

ポテンシャルエネルギー曲面（Potential Energy Surface; PES）上の停留点（stationary points）、すなわち極小点（minima）および鞍点（saddle points）の探索と特性評価は、化学反応動力学や分光学的な分子物性を理解する上で中心的な役割を果たす物理化学的課題である。基底状態における局所的な極小点の決定は、多くの信頼性の高いアルゴリズムの確立により、現代計算化学においては比較的ルーチンな作業と見なされている。しかしながら、鞍点（遷移状態）の探索は、小規模な分子系であっても依然として非自明かつ困難な課題である。特に、鞍点構造においては化学結合が著しく伸長している場合が多く、単一参照（Single Reference; SR）法では、その電子状態を正確に記述するために必要な静的相関（static correlation）を十分に捉えきれない場合がある。

さらに、励起状態のPESに関する調査は、基底状態と比較して量子化学的手法に対する要求が厳しく、極小点の探索でさえも組織的な研究例は限られている。このような臨界的な状況、すなわち結合解離や励起状態においては、多参照（Multireference; MR）法が必須となる。具体的には、多参照配置間相互作用（MR-CISD）、完全活性空間（CAS）、二次摂動論（CASPT2）、そして多参照平均結合対汎関数（MR-ACPF）や多参照平均二次結合クラスター（MR-AQCC）などが挙げられる。特にMR-ACPF/AQCC法は、MR-CI法に対して変分的なアプローチに基づきつつ、サイズ示量性（size-extensivity）の補正を可能にする手法として重要である。

幾何構造探索において、解析的勾配（analytic gradients）の利用可能性は計算効率を左右する決定的な要因である。SR法とは対照的に、MR法における解析的勾配の実装は一般的ではない。そのため、多くの構造最適化は多配置自己無撞着場（MCSCF）レベルに留まることが多く、動的相関（dynamic correlation）の主要部分が欠落することで、構造やエネルギー論に重大なアーティファクト（人為的誤差）が生じるリスクがあった。

本稿で取り上げるDallosらの研究（2002年）の背景には、COLUMBUSプログラムシステムにおいてMR-CISD法の効率的な解析的勾配が開発されたという技術的進展がある。MR-ACPF/AQCC法も変分的な性質を持つため、MR-CISDの形式論と密接な類似性をもって解析的勾配を計算することが可能である。従来、この勾配形式論は単一状態MCSCF（single-state MCSCF）に基づくMR-CISD計算に限定されていたが、Dallosらの研究の直前に、状態平均MCSCF（state-averaged MCSCF）に基づくMR-CISD勾配法へと拡張され、励起状態のバランスの取れた記述が可能となった。

本稿では、Quappらによって開発されたReduced Gradient Following (RGF) 法をCOLUMBUSシステムへ実装し、基底状態および励起状態のPES上での大域的な停留点探索を実現した手法について、その数理的背景から実装の詳細、そして実際の適用例に至るまでを包括的に解説する。

1. Reduced Gradient Following (RGF) 法の数理的定式化#

RGF法は、Distinguished Coordinate Method（区別された座標法）の一般化として位置づけられる手法であり、停留点を特徴づける勾配条件を緩和（reduce）することで、停留点間を連結する曲線を定義し、それを追跡するアルゴリズムである。

1.1 基本方程式と射影演算子#

通常の停留点探索における条件は、エネルギー $E(x)$ の勾配ベクトル $\nabla E(x)$ がゼロになることである。

\nabla E(x) = 0

RGF法において、この条件は一つの方程式を除いて緩和される。Quappらによる本来の形式では、特定の座標 $x^k$ を除く成分について勾配がゼロになる条件が課されたが、これは任意の固定された探索方向 $r$ に対して、以下の射影方程式として一般化される。

P_r \nabla E(x) = 0 \quad \cdots (1)

ここで、 $P_r$ は探索方向 $r$ に直交する部分空間への射影行列（projector matrix）であり、 $P_r r = 0$ を満たす。もし探索方向が特定の座標軸（例えば $k$ 番目の座標）に沿っている場合、この射影行列 $P$ は、 $(N-1) \times N$ の行列として以下のように表現される。

P = \begin{pmatrix} 1 & 0 & \cdots & 0 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 0 & \cdots & 1 \end{pmatrix}

この行列は、単位行列から $k$ 番目の行（探索方向に対応）を取り除いた形をしている。式(1)を満たす点 $x$ の集合は、PES上の停留点同士を結ぶ曲線を定義する。ある既知の停留点（例えば極小点）から出発し、この曲線（RGF曲線）を追跡することで、目的とする鞍点に到達することが可能となる。

1.2 曲線の追跡：Predictor-Corrector法#

RGF曲線を数値的に追跡するために、Predictor-Corrector（予測子-修正子）法が採用されている。

Predictor（予測）ステップ#

RGF曲線上の点 $x(t)$ において、式(1)をパラメータ $t$ で微分することで、曲線の接線ベクトル（tangent vector） $x'$ を求める方程式が導出される。

\frac{d}{dt} P_r \nabla E(x(t)) = P_r H(x(t)) \frac{dx(t)}{dt} = 0

すなわち、

P_r H(x(t)) x' = 0 \quad \cdots (2)

ここで $H(x(t))$ は点 $x(t)$ におけるHessian行列（力の定数行列）である。式(2)は $N-1$ 個の線形方程式系であり、QR分解を用いて解くことができる。これにより得られた接線ベクトル $x'_m$ を用いて、次の点 $x_{m+1}$ を予測する。

x_{m+1} = x_m + \sigma \frac{L}{\|x'_m\|} x'_m \quad \cdots (3)

ここで $L$ はステップ長（steplength）パラメータであり、 $\sigma$ は探索の方向（符号）を制御する。

Corrector（修正）ステップ#

予測された点 $x_{m+1}$ は、一般に厳密には式(1)の条件 $P_r \nabla E(x) = 0$ を満たさない。そのため、Newton-Raphson法に類似した反復法を用いて、射影勾配がゼロになる部分空間へと座標を修正する。この修正ステップを経ることで、計算点は常にRGF曲線近傍に保たれ、安定して鞍点へと誘導される。最終的に鞍点近傍に到達した後は、より収束精度の高いGDIIS（Geometry Direct Inversion in the Iterative Subspace）法に切り替えて構造の精密化を行う。

2. 実装における計算戦略とHessian行列の取り扱い#

多参照法を用いた計算、特にMR-CISDやMR-AQCCレベルでの計算は、単一参照法に比べて計算コストが極めて高い。RGF法の実装において最もクリティカルな問題は、接線ベクトルの算出（式2）およびCorrectorステップにおいてHessian行列 $H$ の情報が必要となる点である。

2.1 MR法におけるHessian計算の困難性#

現状のCOLUMBUSの実装において、MR-CISD/AQCCの解析的二次微分（Hessian）は利用不可能である。これは形式論の複雑さに起因する。したがって、Hessianは解析的勾配の有限差分（finite differences of analytic gradients）として数値的に計算する必要がある。自由度 $N$ の分子系において、完全な数値Hessianを求めるには、対称性を考慮しない場合 $6N$ 回程度の勾配計算が必要となり、これは大域的な探索において計算時間的に許容しがたいコストとなる。

2.2 Hessian更新（Update）手法の導入#

Dallosらは、計算効率と探索の安定性のバランスを最適化するために、Hessianの取り扱いについて以下の3つのケース（Case）を比較検討している。

Case 1: 完全再計算 (Reference)
- RGF探索の各ステップでHessian行列を数値的に再計算する。
- 最も計算コストが高いが、最も正確な情報が得られる。ベンチマークの基準として使用される。
Case 2: 初期計算 + 更新 (Most Efficient)
- 探索の開始点（最初の1点）でのみHessianを計算する。
- 以降のすべてのステップ（PredictorおよびCorrector）では、Bofillによって開発された更新式を用いてHessianを推定する。
- 最も計算コストが低い手法である。
Case 3: Predictorでの再計算 + Correctorでの更新 (Hybrid)
- 幾何構造の変化が大きいPredictorステップにおいてはHessianを再計算する。
- 微修正を行うCorrectorステップでは更新式を用いる。
- Case 1とCase 2の中間的な手法であり、更新のみでは収束しない困難なケース向けである。

2.3 Hessian計算のレベル#

さらに、Hessian計算自体に用いる量子化学計算のレベル（基底関数や理論レベル）を、勾配計算（MR-CISD/AQCC）よりも下げる近似も有効である。RGF法は大域的な探索手法であり、経路の途中で厳密なHessianは必ずしも必要ではないため、Hessian計算にはより小さな基底関数や、SCF/MCSCFレベルの計算を用いることでコストを削減する戦略が取られた。

検証の結果（論文Table 1）、Bofillの更新式を用いた Case 2 が、反復回数を著しく増加させることなく探索を完了できることが示された。例えば、アセチレンの $cis \to trans$ 異性化の探索において、Case 1（完全計算）では25ステップを要したのに対し、Case 2（更新）では18ステップで収束するなど、更新式がHessianの局所的な曲率変化を適切に捉えていることが実証された。

3. 適用事例研究1：ホルムアルデヒド ( $H_2CO$ ) の $S_0$ 曲面#

最初の適用例として、ホルムアルデヒドの基底状態（ $S_0$ ）における停留点探索が行われた。これには、解離反応 ( $H_2 + CO$ ) およびヒドロキシカルベン ( $HCOH$ ) への異性化反応が含まれる。

3.1 計算詳細#

手法: MR-CISD および MR-AQCC
基底関数: cc-pVTZ
活性空間:
- MCSCF計算にはREDVAL（Reduced Valence）空間を使用。
- MR-CISD/AQCC計算には、計算コスト削減のためREDVAL-2ex空間（主要な4軌道からの単・二電子励起に限定した参照空間）を採用。
- 内殻軌道（ $1a_1, 2a_1$ ）は凍結（frozen core）。

3.2 探索された反応経路と構造#

図1（論文 Figure 1）に示されるように、大域的な極小点であるホルムアルデヒド構造 ( $M_1, C_{2v}$ ) を起点として、以下の反応経路が探索された。

分子解離: $H_2CO \to H_2 + CO$
- 遷移状態 $F_4 (C_s)$ を経由。
- バリア高さはMR-AQCCで87.1 kcal/molと算出され、FellerらによるCCSD(T)の推定値（87.4 kcal/mol）と極めて良い一致を示した。
- 生成系 ( $M_2$ ) への反応エネルギーは4.5 kcal/molと算出された。
ヒドロキシカルベンへの異性化: $H_2CO \to t-HCOH (M_3) / c-HCOH (M_4)$
- トランス体 ( $M_3$ ) へは遷移状態 $F_2$ を経由し、シス体 ( $M_4$ ) へは遷移状態 $F_1$ を経由する。
- $t-HCOH$ ( $M_3$ ) の相対エネルギーは52.4 kcal/mol、 $c-HCOH$ ( $M_4$ ) は58.0 kcal/molであった。

3.3 鞍点 $F_2$ の次数に関する議論#

特筆すべき成果として、遷移状態 $F_2$ （ホルムアルデヒド $\to$ トランス-ヒドロキシカルベン）の特性評価が挙げられる。過去のQuappらによるSTO-3G/SCFレベルの研究では、 $F_2$ は二次（second-order）の鞍点であると報告されていた。しかし、本研究におけるMR-CISD/AQCCレベルの計算では、 $F_2$ は明確に一次（first-order）の鞍点であることが確認された。これは、鞍点の特性（虚振動数の数）が計算レベルに強く依存する場合があることを示す重要な事例であり、高精度なMR法による検証の必要性を裏付けている。

3.4 幾何構造の比較#

表2および表4に示されるように、REDVAL-2ex-AQCCレベルで最適化された幾何構造は、既存のCCSD(T)や実験値と良好に一致している。例えば、 $F_4$ 遷移状態のC-O結合長は1.171 Å（MR-AQCC）であり、CCSD(T)の1.163 Åと非常に近い。また、ヒドロキシカルベン異性体の構造についても、以前のDFTやQCISD計算と比較して妥当な結果が得られており、MR法が単一参照法と同等以上の精度で構造決定を行えることが示された。

4. 適用事例研究2：アセチレン ( $C_2H_2$ ) の励起状態 ( $T_1, T_2$ )#

二つ目の適用例として、アセチレンの低次三重項状態（ $T_1, T_2$ ）におけるPESの探索が行われた。これらの状態は平面構造および非平面構造を含み、複雑なトポロジーを持つため、多参照性が強く現れる系である。

4.1 計算詳細と状態平均#

手法: MR-AQCC（MCSCF軌道に基づく）
基底関数: cc-pVTZ
状態平均（State-Averaging）:
- $S_0, T_1, T_2$ を含む複数の電子状態をバランスよく記述するため、状態平均MCSCF計算（ $S_0$ + 3つの三重項 + 2つの一重項）を実行し、その軌道を用いてMR計算を行った。
- 参照空間はCAS(4,4)（ $1\pi_u, 1\pi_g$ 軌道に4電子）から構築。

4.2 探索結果：平面および非平面鞍点#

アセチレンの $T_1$ ( $1^3B_2$ ) および $T_2$ ( $1^3A_u$ ) 状態において、cisおよびtrans型の極小構造が存在する。本研究では、これらを連結する鞍点として、以下の二つのタイプが特定された。

平面鞍点 (Planar Saddle Points)
- $C_s$ 対称性を持つ構造。
- $T_1$ 状態における平面鞍点は二次の鞍点（second-order）であることが判明した。
- $T_2$ 状態における平面鞍点 ( $1^3A''$ ) は一次の鞍点であり、cis構造に対するバリアは0.30 eVと算出された。
非平面鞍点 (Out-of-Plane Saddle Points)
- $C_2$ 対称性を持つ構造。
- $T_1$ 状態 ( $1^3B$ ) において、バリア高さ0.62 eVの一次鞍点が特定された。
- $T_2$ 状態 ( $1^3A$ ) においても、バリア高さ0.26 eVの一次鞍点が存在する。

4.3 多参照法の必要性と既存研究との比較#

アセチレンの励起状態においては、SR法（例えばSR-CISDやEOM-CCSD）を用いた場合、波動関数の不安定性や記述の破綻が生じやすいことが知られている。特にCuiらによる先行研究では、CASSCFやEOM-CCSDを用いた非平面構造の研究が行われているが、CASSCFでは動的相関が欠落し、EOM-CCSDでは参照波動関数の単一参照性が制限となる。本研究の結果（表5）は、 $T_1$ 状態の非平面鞍点において、SR-CISDの結果（Vacek et al.）と比較して結合長が長くなる傾向（MR効果による結合のより適切な記述）を示しており、MR-AQCC法による記述がより物理的に妥当であることを示唆している。また、MR-AQCCにより算出された cis- $T_1$ 構造の励起エネルギー（3.93 eV）は、SherrillらによるCCSD(T)の推定値（3.91 eV）と極めて良く一致し、実験的な推定値（3.58 eV）が低すぎるという彼らの主張を支持する結果となった。

特に、 $T_2$ 状態の非平面鞍点 ( $1^3A$ ) は、さらに高次の三重項状態 ( $2^3B$ ) との交差（crossing）に関与しており、アセチレンの光解離ダイナミクスを理解する上で重要である。このような縮退付近の領域では多参照法によるアプローチが不可欠であり、本研究はMR法による最初の完全な構造最適化の例として位置づけられる。

5. 結論と総括#

本論文は、RGF法を多参照量子化学計算プログラムCOLUMBUSに実装することで、基底状態のみならず励起状態のポテンシャルエネルギー曲面における大域的な停留点探索が可能になったことを報告している。

手法の確立: 射影勾配方程式 $P_r \nabla E = 0$ に基づくRGF法は、MR-CISD/AQCCレベルの計算においても、BofillのHessian更新法を併用することで、計算コストを抑制しつつ安定して鞍点を特定できることが実証された。
精度の向上: ホルムアルデヒドの系において、従来の手法（STO-3G/SCF等）で見られた鞍点の次数に関する誤り（ $F_2$ 構造）を正し、高精度な幾何構造パラメータとエネルギー障壁を提供した。
励起状態への展開: アセチレンの $T_1, T_2$ 状態において、平面および非平面の鞍点を網羅的に探索し、多参照性が支配的な領域における構造決定において本手法が強力なツールとなることを示した。

MR法における解析的二次微分の欠如という制約を、Hessian更新アルゴリズムとRGF法の組み合わせによって克服した本研究は、光化学反応や複雑な電子状態を持つ系の反応経路探索における重要なマイルストーンであると言える。

参考文献#

M. Dallos, H. Lischka, E. Ventura do Monte, M. Hirsch, W. Quapp, “Determination of Energy Minima and Saddle Points Using Multireference Configuration Interaction Methods in Combination with Reduced Gradient Following: The $S_0$ Surface of $H_2CO$ and the $T_1$ and $T_2$ Surfaces of Acetylene”, Journal of Computational Chemistry, 23, 576-583 (2002).
W. Quapp, M. Hirsch, O. Imig, D. Heidrich, Journal of Computational Chemistry, 19, 1087 (1998).
H. Lischka, M. Dallos, R. Shepard, Molecular Physics, 100, 1647 (2002).
J. M. Bofill, Journal of Computational Chemistry, 15, 1 (1994).
D. Feller, M. Dupuis, B. C. Garrett, Journal of Chemical Physics, 113, 218 (2000).
Q. Cui, K. Morokuma, J. F. Stanton, Chemical Physics Letters, 263, 46 (1996).