幾何構造最適化における初期ヘシアン推定の数理的基礎と拡張：Schlegel Hessianについて

4023 words

20 minutes

幾何構造最適化における初期ヘシアン推定の数理的基礎と拡張：Schlegel Hessianについて

2026-01-25

Computational Chemistry

Geometry Optimization

Hessian

Force Constants

Algorithm

Schlegel

last_modified: 2026-01-25

生成AIによる自動生成記事に関する免責事項: 本記事は、H. Bernhard Schlegelらによる原著論文 “Estimating the Hessian for gradient-type geometry optimizations” (1984) および “Estimating stretching force constants for geometry optimization” (1997) の内容に基づき、大規模言語モデルによって作成された解説記事です。記述内容は提示された文献の範囲内において正確性を期していますが、厳密な証明や詳細な実装については、原著論文および関連する計算化学の専門書を参照してください。

1. 序論：勾配法による構造最適化とヘシアンの役割#

計算化学における幾何構造最適化は、ポテンシャルエネルギー曲面（PES）上の停留点探索問題であり、現代の多くの量子化学計算プログラムにおいて中心的なタスクの一つである。この最適化プロセスの効率、すなわち平衡構造への収束速度は、初期構造の質とともに、エネルギーの二次微分行列であるヘシアン（Hessian、力定数行列）の初期推定精度に強く依存する。

ヘシアンの初期推定に関して、両極端な2つのアプローチが存在する。一つは、初期ヘシアンとして単位行列（Unit Matrix）を用いる方法である。これは偏りのない選択ではあるものの、分子の結合性や原子の種類といった有用な構造情報をすべて破棄してしまうため、効率的ではない。特に、原子価座標（結合伸縮や結合角など）間の結合（カップリング）が無視されるため、環状分子のように座標同士が本質的に強く結合している系においては、最適化の収束が著しく損なわれる。

もう一つの極端なアプローチは、初期段階で解析的なヘシアンを厳密に計算することである。ポテンシャル面が二次形式で近似できる場合、正確なヘシアンを用いればわずか1ステップで最適解に到達可能である。しかし、 $N$ 変数系における解析的ヘシアンの計算は、勾配計算の約 $N$ 倍の計算コストを要するため、大規模な系においては計算資源的に正当化できない場合が多い。

したがって、これらの極端なアプローチに対する折衷案としての計算コストを抑えつつ、分子の化学的性質（原子価力場）を反映した推定手法の開発に需要があると見込まれた。H. Bernhard Schlegelは1984年、単純な原子価力場（Valence Force Field）を用いてヘシアンを経験的に構築し、それを最適化空間へ変換する手法（通称Schlegelヘシアン）を提唱した。本稿では、この手法の理論的枠組みと、1997年に行われた重元素へのパラメータ拡張について詳述する。

2. 理論的枠組み：冗長原子価座標によるヘシアン構築#

Schlegelの方法論の中核は、最適化に使用する座標系（非冗長な内部座標やZ-matrix）に依存せず、分子のトポロジーに基づいて定義される「冗長原子価座標（Redundant Valence Coordinates）」上で初期力定数を推定する点にある。

（注：Schlegelの方法論の中核は、最適化に用いられる内部座標系の選択に依存せず、分子のトポロジーに基づいて初期力定数を生成できる点にある。ただし、最終的な最適化挙動や収束効率は、選択された内部座標系の性質に依存する。）

2.1 冗長原子価座標系の定義#

まず、原子のデカルト座標から、結合距離、結合角、二面角（捩れ角）、平面外変角（Out-of-plane bend）などの原子価座標を生成する。原子間の結合の有無は、原子間距離と共有結合半径（covalent radii）に基づいて判定される。具体的には、原子間距離が2つの原子の共有結合半径の和の1.35倍未満である場合、それらの原子は結合しているとみなされる。

この座標系は、分子の自由度（ $3N-6$ ）よりも多くの変数を持つため「冗長（redundant）」と呼ばれるが、分子の局所的な化学的性質（結合の強さや硬さ）を記述するには最適である。

2.2 力定数の経験的推定式#

原子価座標系における対角力定数 $F$ は、分光学における経験則を計算化学用に簡略化した式を用いて推定される。

2.2.1 結合伸縮（Bond Stretch）：Badger則の適用#

結合伸縮の力定数 $F_{str}$ は、Badger則として知られる以下の形式で近似される。

F_{str} = \frac{A}{(r - B)^3} \quad (1)

ここで、 $r$ は結合距離、 $A$ と $B$ は定数である。定数 $A$ はすべての結合に対して一定の値が用いられる。一方、定数 $B$ は結合する2つの原子が周期表のどの「周期（Row）」に属しているかにのみ依存する。この式により、結合距離が短い（強い結合）ほど力定数が大きくなるという物理的直観が反映される。

2.2.2 結合角（Angle Bend）#

結合角の力定数 $F_{bend}$ については、伸縮ほど一般的な法則が存在しないため、1984年の段階では簡便な定数割り当てが採用された。

重原子のみで構成される結合角： $A = 0.250$ (hartree/rad $^2$ )
水素原子を含む結合角： $A = 0.160$ (hartree/rad $^2$ )

これにより、水素を含む角度変形の方が柔らかいという傾向が表現される。

2.2.3 捩れ角（Torsion）#

捩れ（二面角）の力定数 $F_{tors}$ は、結合の多重度を考慮するために、結合軸となる原子間の距離 $r$ に依存する形式が用いられる。二重結合周りの回転障壁は単結合周りよりも一桁大きくなる可能性があるためである。

F_{tors} = A - B(r - r_{cov}) \quad (r < r_{cov} + A/B) \quad (2)

ここで、 $r_{cov}$ は共有結合半径の和（単結合距離の参照値）である。この式は、結合距離が短くなる（二重結合性が増す）につれて力定数が増大することを表現している。具体的なパラメータとして、 $A=0.0023$ , $B=0.07$ が用いられる。

2.2.4 平面外変角（Out-of-plane bend）#

平面性を持つ中心原子（例： $sp^2$ 炭素やカルボニル基）からのピラミッド化に対する抵抗力を記述するために、平面外変角座標が導入される。その力定数 $F_{oop}$ は、構造の非平面度 $d$ に依存して減衰するように設計されている。

F_{oop} = A d^4 \quad (3)

d = 1 - \frac{\mathbf{r}_1 \cdot (\mathbf{r}_2 \times \mathbf{r}_3)}{|\mathbf{r}_1| |\mathbf{r}_2| |\mathbf{r}_3|} \quad (4)

ここで $\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3$ は中心原子から伸びる3つの結合ベクトルである。平面構造においては $d=1$ となり最大の力定数を持つが、構造がピラミッド化する（非平面化する）につれて力定数は急速に減衰する。(式4は一例であり、実装の仕方によって変わることに注意されたし。 $d$ は中心原子の非平面性を表す幾何学的指標であり、実装では結合ベクトルから構成されるスカラー量として評価される。)

3. 座標変換のアルゴリズム：原子価空間から最適化空間へ#

上述の原子価座標系で推定された力定数行列 $F^{val}$ は、通常対角行列（あるいは近接相互作用のみを含む疎行列）であるが、これを実際の最適化計算で使用される座標系へ変換する必要がある。この変換は、一度デカルト座標系を経由して行われる。

3.1 原子価座標からデカルト座標への変換#

原子価座標の微小変位 $q'$ とデカルト変位 $x$ の関係は、Wilsonの $B$ 行列（ここでは $B'$ と表記）によって記述される。

q' = B' x \quad (5)

これを用い、原子価座標上の力定数行列 $F^{val}$ は、以下の変換式によってデカルト座標上の力定数行列 $F^{cart}$ へ変換される。

F^{cart} = B' F^{val} {B'}^T \quad (6)

この操作により、原子価空間での局所的な剛性情報が、原子間の空間的な接続情報としてデカルト座標上に展開される。 $F^{val}$ が対角行列であれば、 $B'$ 行列を明示的に保存する必要がなく、計算は特に単純化される。

3.2 デカルト座標から非冗長内部座標への変換#

最適化が、Z-matrixなどで定義された非冗長な内部座標 $q$ で行われる場合、さらに変換が必要となる。非冗長座標とデカルト座標の関係を $q = B x$ とする（ここでの $B$ は非冗長系に対する行列）。

力定数の変換には $B$ の逆行列が必要となるが、 $B$ は正方行列ではない（デカルト座標の自由度 $3N$ に対し、内部座標は $3N-6$ ）ため、一般化逆行列が用いられる。

B^{-1} = (B^T M B)^{-1} B^T M \quad (7)

ここで $M$ は補助行列であり、並進・回転の自由度を処理するための重み付け行列である。

最終的な内部座標系での力定数行列 $F^{int}$ は以下のように求められる。

F^{int} = B^{-1} F^{cart} (B^{-1})^T + \sum_k \frac{\partial (B^{-1})_{ik}}{\partial q_j} f^{cart}_k \quad (8)

第2項は座標系の曲線性（curvilinear nature）に由来する項(共変微分と等価と考えられる。)であるが、初期ヘシアンの推定においては $F^{cart}$ 自体が近似値であるため、この項は通常無視される。したがって、実用的な変換式は以下の通りとなる。

F^{int} \approx B^{-1} F^{cart} (B^{-1})^T \quad (9)

この一連の変換プロセス（Valence $\to$ Cartesian $\to$ Internal）を経ることで、ユーザーが定義した最適化座標系の良し悪しに関わらず、分子固有の物理的接続情報を反映したヘシアンを得ることが可能となる。

4. パラメータの拡張と精密化（1997年の展開）#

1984年の提唱以降、計算化学の適用範囲は飛躍的に拡大し、第3周期以降の重元素を含む計算が日常的に行われるようになった。これを受け、1997年にWittbrodtとSchlegelは、第6周期（Cs-At）までの元素に対応するためのBadger則パラメータの再決定を行った。

4.1 データセットと計算手法#

パラメータ決定のために、第1周期から第6周期までの元素を含む80以上の単純な分子（二原子分子、水素化物、ハロゲン化物など）が選定された。これらの分子に対し、以下のレベルで構造最適化と周波数解析（解析的力定数の算出）が行われた。

HF/3-21G: 第1～5周期
HF/LanL2DZ: 第6周期を含む系
MP2/6-31G*: 第1～3周期
B3LYP/6-31G*: 第1～3周期およびBr
B3LYP/LanL2DZ: 第4周期以降、および第5, 6周期

4.2 Badger則パラメータ $B$ の再評価#

1984年の研究との互換性を保つため、定数 $A$ は $1.734$ (hartree/bohr) に固定された。その上で、計算された結合伸縮力定数 $F_{str}$ と結合距離 $r$ から、各周期ペアに対応するパラメータ $B$ が逆算された。

B = r - \left( \frac{1.734}{F_{str}} \right)^{1/3} \quad (10)

得られた結果（Table 1, 2）から、パラメータ $B$ は結合する原子の周期に対して対数的な依存性を示す傾向があることが判明した。ただし、アルカリ金属水素化物（LiH, NaH等）においては、例外的に線形な依存性が観測された。

4.3 精度検証と例外#

再パラメータ化されたBadger則による推定力定数は、B3LYPレベルで計算された真の力定数に対し、相関係数0.964という良好な一致を示した。平均絶対相対誤差は11.9%であり、最適化の初期推定として要求される精度（10-20%）を十分に満たしていることが確認された。

しかし、いくつかの重要な例外も報告されている。

$Ni(CO)_4$ : 推定誤差が40%と最大であった。これは配位結合（dative bonding）の性質が強く、Badger則の単純なモデルでは記述が困難であるためである。パラメータ $B$ の値（1.199 bohr）も特異的に小さかった。
アルカリ金属水素化物: LiHやNaHなどは $B$ パラメータが極端に小さく、全体的な平均化からは除外された。
非共有結合相互作用: 本手法は水素結合やファンデルワールス相互作用については考慮していない。

5. 議論：座標系の選択と最適化の効率#

Schlegelヘシアンの利点は、適切な初期値を与えることにあるが、最適化の収束は座標系の選択にも依存する。1984年の論文では、座標間のカップリング（Coupling）を減らすことの重要性が議論されている。

特に問題となるのは、剛直なモード（stiff mode、結合伸縮など）と柔軟なモード（flexible mode、結合角変形など）が強く結合してしまうケースである。例えば、6員環の平面構造を記述する際、ある角度の変化が環を閉じる結合の距離に直接影響を与えるような座標定義を行うと、ヘシアンの非対角項が大きくなり、最適化はポテンシャル面上の「長く狭い谷」を進むような困難な挙動を示す。 Schlegelヘシアンを用いることで、こうした座標間の物理的なカップリングを初期ヘシアンに埋め込むことができるため、環状分子などの「病的な」系に対しても、単位行列を用いる場合と比較して格段に優れた収束性を実現できる。

6. 結論#

H. Bernhard Schlegelによって確立され、その後Wittbrodtらによって拡張されたヘシアン推定法は、以下によって特徴づけられる。

物理的基盤: 原子の接続性と周期表に基づく単純な経験則（Badger則）を利用し、計算コストをかけずに妥当な初期力場を構築する。
アルゴリズムの頑健性: 冗長原子価座標を経由する変換アルゴリズムにより、どのような最適化座標系（Z-matrix等）を選択しても、分子のトポロジーを反映したヘシアンを提供できる。
拡張性: 第6周期までの重元素を含む広範な化学種に対し、10-20%程度の誤差範囲で力定数を推定可能である。

この手法は、解析的ヘシアン計算のコストを回避しつつ、勾配法による構造最適化のステップ数を削減する実利的な解法として、今日でも広く利用されている。

参考文献#

Schlegel, H. B. (1984). Estimating the Hessian for gradient-type geometry optimizations. Theoretica Chimica Acta (Berl.), 66, 333-340.
Wittbrodt, J. M., & Schlegel, H. B. (1997). Estimating stretching force constants for geometry optimization. Journal of Molecular Structure (Theochem), 398-399, 55-61.