Hessian・振動数・熱化学

Hessian・振動数・熱化学#

勾配は平衡にあるかどうかを教えました。Hessian — エネルギーの 2 次微分の行列 — は、それがどの種類の停留点かを教え、化学者が実際に極小で測る 2 つのもの — 振動数（赤外/ラマンスペクトル）と熱化学（平衡と速度を予測する自由エネルギー）— を解き放ちます。これは実在し検証済みの qc-rs 機能（m.scf.hessian, qc.thermo.frequencies）で、本マニュアルの他の場所ではまだ扱っていません。

理論:なぜ 2 次微分に軌道応答が要るのか#

もう一段エネルギーを微分すると、何かが変わります。勾配 — 1 次微分 — は軌道応答項を必要としませんでした。 Hellmann–Feynman がそれを消し去るからです:SCF は軌道について停留している（\(\mathbf g=0\)、SCF 収束理論）ので、 \(\partial E/\partial R_A\) は演算子の明示的な \(R\) 依存性しか見ません。その停留条件をもう一度微分すると、軌道自体の応答が再び現れます — なぜなら、ある核が動くときに軌道がどう緩和するかは、他のすべての核の位置の関数だからです:

\[ \frac{\partial^2 E}{\partial R_A\,\partial R_B} = \underbrace{E^{[2]}_{AB}}_{\text{骨格:2 次微分積分}\;\times\;D,\,W} \;+\; \underbrace{\sum_{ai} \frac{\partial F_{ai}}{\partial R_A}\,U^{B}_{ai} \;+\; \text{c.c.}}_{\text{軌道応答 — }U^{B}\text{が必要}}. \]

骨格項は易しい半分です — 2 次微分積分（\(\partial^2 T\), \(\partial^2 V_{ne}\), \(\partial^2(\mu\nu|\lambda\sigma)\), \(W\) 経由の \(\partial^2 S\)）を無摂動密度と縮約するだけで、勾配とまったく同じですが微分次数が 1 つ上です。軌道応答項は新しい物理です:\(U^B = \partial\mathbf C/\partial R_B\)、核 \(B\) が動くとき軌道自体がどう緩和するか、**結合摂動 Hartree–Fock/Kohn–Sham（CPHF/CPKS）**方程式

\[ (\mathbf A+\mathbf B)\,U^B = -\,\mathbf b^B, \]

を解いて求めます。これは SCF 収束理論の拡大 Hessian 機構と同じ軌道 Hessian 演算子 \(\mathbf A+\mathbf B\) です — ここでは固有値問題ではなく線形連立方程式として再利用され、右辺 \(\mathbf b^B\) は 1 次微分 Fock \(F^{B,x}\) から幾何微分として作られます。これが勾配を超える唯一の真に新しい要素です:qc-rs の CPHF/CPKS 応答エンジンは既に存在する（安定性解析も駆動する）ので、Hessian はその上に「骨格 + fold-back」のレシピを組み立てます:

        flowchart LR
    S["骨格項<br/>2 次微分積分 × D, W"] --> SUM["+"]
    G["1 次微分 Fock F^B,x<br/>（勾配用に既に構築済み）"] --> RHS["CPHF 右辺 b^B"]
    RHS --> CPHF["(A+B) U^B = −b^B を解く<br/>（SCF 応答エンジンを再利用）"]
    CPHF --> FOLD["Fold-back 項<br/>Σ ∂F/∂R_A · U^B"]
    FOLD --> SUM
    SUM --> H["分子 Hessian ∂²E/∂R_A∂R_B"]

アルゴリズム 2 (分子 Hessian の組み立て)

入力: 収束した SCF（\(\mathbf C\), \(\varepsilon\), \(\mathbf D\)）、CPHF 応答エンジン。 出力: Hessian \(\mathbf H \in \mathbb R^{3n_{\text{atom}}\times 3n_{\text{atom}}}\)。

骨格ブロック — 一電子、二電子（4 中心または RI）、（KS なら）XC の 2 次微分項 — を、固定密度 \(\mathbf D\) とエネルギー重み付き密度 \(\mathbf W\) と縮約して積み上げる。
幾何微分 Fock \(F^{B,x}\) と重なり \(S^{B,x}\) をすべての核 \(B\) について構築する（勾配が既に必要とする 1 次微分積分を再利用）。
各摂動 \(B\) について、\((\mathbf A+\mathbf B)U^B = -\mathbf b^B\) を解いて軌道応答 \(U^B\) を求める（CPHF/CPKS の求解 — 核ごとに 1 回の線形解、合計 \(3n_{\text{atom}}\) 回）。
\(U^B\) を fold back する:\(\sum_{ai}(\partial F_{ai}/\partial R_A)\,U^B_{ai} + \text{c.c.}\) を骨格に加える。
古典的な核間反発 Hessian \(\partial^2 V_{nn}/\partial R_A\partial R_B\)（閉形式）を加える。
自己検査:音響和則 \(\sum_B \mathbf H_{AB} \approx \mathbf 0\) が成り立たねばならない（分子全体の剛体並進はどの 2 次微分も変えない）— 勾配の並進不変性チェックの Hessian 版。

Hessian からスペクトルへ:基準振動#

デカルト座標の Hessian は並進・回転・真の振動を混ぜています。それらを解きほぐすのは、短い標準的な手順です:

Hessian を質量重み付けする、\(\tilde H_{Ai,Bj} = H_{Ai,Bj}/\sqrt{m_A m_B}\)（質量は原子単位）。
剛体並進・回転を射影で除く — 非線形分子で 6 方向、線形分子で 5 方向（分子軸周りの回転はエネルギーを持たない）— 質量重み付き並進/回転射影子を使う。
射影した \(\tilde H\) を対角化する:各固有値 \(k\) が基準振動の力の定数で、対応する固有ベクトルが 基準振動（その振動の原子変位パターン）。
振動数へ変換:\(\tilde\nu\,[\text{cm}^{-1}] = \sqrt{k}\cdot\text{（単位換算）}\)。負の \(k\) は虚数振動数を与えます — 真の極小ではなく鞍点の署名です（反応座標に沿った 1 つの虚数モードは、まさに遷移状態の見た目です）。

だから \(N\) 原子の非線形分子は \(3N-6\) 個の振動モードを持ち、線形分子は \(3N-5\) 個（除く回転が 1 つ少ない）です。水（\(N=3\)、非線形）は 3 個;H₂（\(N=2\)、線形）は 1 個 — どちらも下で検証します。

熱化学:振動数から自由エネルギーへ#

調和振動数は標準的な理想気体統計力学（並進 + 回転 + 振動の分配関数）に入力され、化学者が実際に欲しい量を、与えた温度・圧力で生成します:

零点エネルギー \(\text{ZPE} = \tfrac12\sum_k h c\,\tilde\nu_k\) — \(T=0\) でも各調和モードが寄与する振動基底状態エネルギー。
エンタルピー \(H\)、エントロピー \(S\)、Gibbs 自由エネルギー \(G=H-TS\)、熱容量 \(C_v,C_p\) — 各々、並進・回転・振動の分配関数から通常どおり構築され、電子エネルギーに加えられます。

\(G\) こそが化学平衡を実際に予測し、遷移状態理論を通じて反応速度も予測します — 調和振動数は「1 つの構造での 1 つの SCF エネルギー」からそこへの橋です。

使い方#

Hessian は、勾配と同じく SCF ステップの結果です:

import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r").run()

H = np.asarray(m.scf.hessian)      # 形 (3*natom, 3*natom) = 水なら (9, 9)
H.sum(axis=1)                       # どの行も ~0:音響和則の自己検査

qc.thermo.frequencies(mychk) が、質量重み付け・射影・対角化・熱化学まで、1 回の呼び出しで残りを行います:

fr = qc.thermo.frequencies(m)

fr.frequencies     # 調和波数の配列、cm⁻¹
fr.n_imaginary     # 虚数（負の k）モードの数:真の極小では 0
fr.norm_mode       # 基準振動の変位ベクトル
fr.thermo          # dict: temperature, pressure, zpe, e_tot, h_tot, g_tot, s_tot, cv_tot, cp_tot

収束した構造は、実務上だけでなく理論上も前提条件です:真の極小から外れた（\(\mathbf g\ne 0\)）構造で計算した振動数は、意味のある調和振動数ではありません — 実際の解析の前には必ず .opt() （勾配と構造最適化）を走らせてください。

総合例#

import qc, numpy as np
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r").run()

fr = qc.thermo.frequencies(m)
print("frequencies (cm⁻¹):", np.round(np.sort(fr.frequencies), 2))
print("imaginary modes    :", fr.n_imaginary)
print("ZPE (Ha)           :", round(fr.thermo["zpe"], 6))
print("G at 298.15 K (Ha) :", round(fr.thermo["g_tot"], 6))

# frequencies (cm⁻¹): [2041.42 4494.04 4796.72]
# imaginary modes    : 0
# ZPE (Ha)           : 0.025817
# G at 298.15 K (Ha) : 0.033348

# 線形分子:3N-5 = 1 モード、3N-6 ではない
h2 = qc.chk.new(atom="H 0 0 0; H 0 0 0.74", ao="sto-3g", unit="angstrom").scf(ref="r").run()
print("H2 modes:", len(qc.thermo.frequencies(h2).frequencies))   # 1

実数で正の 3 つの振動数と虚数モード 0 個が、水の構造が真の極小であり鞍点でないことを確認します — そしてその数そのもの（非線形の水で 3、線形の H₂ で 1）が、並進/回転射影がうまくいったことの直接的な検査です。

練習 7

同僚が部分的にしか最適化していない構造で振動数を計算しました。n_imaginary が 1 と返ってきます。これは何を意味し、振動数を信じる前に何を確認すべきですか。
H₂（線形、\(N=2\)）はちょうど \(3N-5=1\) 個の振動モードを持つのに、水のような屈曲三原子分子（\(N=3\)、非線形）は \(3N-6=3\) 個です — 同じ原子数の線形分子と比べて、水の「失われた」回転モードはどこに行きましたか。
Hessian は勾配が必要としない CPHF の求解を必要とします。1 次微分がそれを避けられ、2 次微分が避けられないのはなぜか、一文で述べなさい。

解答練習 7

虚数振動数が 1 つあるということは、Hessian が負の固有値を持つ — 構造は極小ではなく鞍点（例:遷移状態）で、不完全な最適化と整合します。どの振動数も本物の振動として信じる前に、.opt() を完全収束まで再実行すべきです（opt.converged と勾配の軌跡を確認）。
線形分子は独立な回転自由度を2 つしか持ちません（自身の軸周りの回転は、その軸周りの慣性モーメントがゼロなのでエネルギーを持たない）ので、6 つの剛体方向のうち 5 つしか除かれず、\(3N-5\) 個の振動が残ります。非線形分子は 3 つの回転すべてを持つので 6 つ除かれ、\(3N-6\) 個が残ります。水は屈曲しているので、除くべき回転自由度が（2 つでなく）完全な 3 つあります;線形三原子分子が水に対して「余分に」持つはずのモードは、まさに水では既にその 3 振動の中に数えられている屈曲運動であり、失われたモードではありません。
Hellmann–Feynman が1 次微分を軌道について停留にするので軌道応答が打ち消し合います;その停留条件をもう一度微分することが、まさに第 2 の核変位に対する軌道自体の応答を再導入することであり、それを与えられるのは CPHF の求解だけです。

エネルギー・力、そして今や振動/熱化学解析が揃ったので、残りのガイド章は環境 — 溶媒和と分散 — へ、そして可視化・ログ・物性スイート全体へ進みます。