# 解析的 Hessian・振動数・熱化学

[勾配の章](analytic-derivatives.md) は、1 次エネルギー微分が軌道応答をまったく必要とせず、収束済み
密度と $W$ だけで済むことを示しました。2 次微分はその単純化を失います — それは
[CPHF/CPKS 機構](linear-response-cphf.md) 経由で計算される真の結合摂動軌道応答を必要とします —
そしてそれが手に入れば、短く標準的なレシピがそれを、化学者が停留点で実際に測定する 2 つのもの:
振動数と熱化学に変換します。本章は完全な骨格+折り返し Hessian 組み立て、真の振動だけを分離する
質量重み付き並進/回転射影、そして qc-rs が評価する厳密な理想気体分配関数の式を導出します —
`.design/76analytic-hessian-frequencies.md` と `python/qc/thermo.py`(PySCF の
`hessian/thermo.py` の直接移植)に基づきます。[ガイドの章](../20-guide/hessian-frequencies-thermo.md)
がコンパクトな使い方ダイジェストを持ちます;本章は完全な導出を持ちます。

## なぜ 2 次微分が軌道応答を再導入するのか

停留条件 $\partial E/\partial C=0$ を**もう一度**微分することが、まさに軌道応答が再登場する場所です
— ある 1 つの原子核が動くにつれて軌道がどう緩和するかは、それ自体が他のあらゆる原子核がどこにある
かの関数であり、その依存性は Hellmann-Feynman がそこでそれを消し去ったがゆえに 1 次微分には見え
なかったのです。完全な 2 次微分はきれいに 2 つの部分に分かれます:

$$
\frac{\partial^2E}{\partial R_A\partial R_B} = \underbrace{E^{[2]}_{AB}}_{\text{骨格(固定 }C\text{)}}
+ \underbrace{\sum_{ai}\frac{\partial F_{ai}}{\partial R_A}\,U^B_{ai} + \text{c.c.}}_{\text{軌道応答 —
}U^B\text{が必要}},
$$

$U^B=\partial\mathbf C/\partial R_B$ は結合摂動方程式 $(\mathbf A+\mathbf B)U^B=-\mathbf b^B$ を解く
ことで求まります — [線形応答の章](linear-response-cphf.md) の $A+B$ 演算子そのものであり、ここでは
固有値問題ではなく線形解法として再利用されます(核 Cartesian 座標ごとに 1 つの右辺、合計
$3n_{\text{atom}}$ 回の解法)。

## 骨格項の完全な形

骨格は「簡単な半分」です — すべての項が*摂動されていない*密度とエネルギー重み付き密度と縮約され、
勾配とまったく同じですが 1 段階高い微分次数であり、勾配の単一の `int2e_ip1` ではなく**3 つ**の
異なる 2 次微分 2 電子積分クラスを使います:

$$
E^{J,[2]}_{AB} - c_xE^{K,[2]}_{AB} = \tfrac12\sum_{\mu\nu\lambda\sigma}\Bigl(D_{\mu\nu}D_{\lambda\sigma} -
\tfrac{c_x}{2}D_{\mu\lambda}D_{\nu\sigma}\Bigr)\,\partial^2_{AB}(\mu\nu|\lambda\sigma),
$$

2 次微分 $\partial^2_{AB}$ は、どの原子対を微分するかに応じて 3 つのブロックパターンに分かれます:

| ブロック | 積分 | 何を微分するか |
|---|---|---|
| 対角(`A,A`) | `int2e_ipip1` | 同じ bra 中心での両微分 |
| 非対角(`A,B`) | `int2e_ip1ip2` | bra 中心で 1 回、ket 中心で 1 回 |
| 対内(`A,B`) | `int2e_ipvip1` | 両方とも bra 対だが異なる bra 殻 |

これに加えて 1 電子/Pulay 骨格(2 次微分の運動エネルギー/核引力積分を $D$ と縮約、2 次微分の重なりを
$W$ と縮約 — 構造的には勾配の Pulay 項の 1 段上)、そして KS 参照では [DFT グリッドの章で導出](dft-xc-quadrature.md)
された同じグリッド機構から構築される XC 2 次微分寄与です。

## CPHF の右辺と折り返し

各原子核の結合摂動解法を駆動する右辺は、*勾配自身の*1 次微分 2 電子積分(`int2e_ip1`、
[既に勾配のために構築済み](analytic-derivatives.md))を再利用します — この部分のために新しく計算
すべきものはなく、再結合するだけです:

$$
F^{A,x}_{\mu\nu} = h^{A,x}_{\mu\nu} + \sum_{\lambda\sigma}D_{\lambda\sigma}\Bigl[(\mu\nu|\lambda\sigma)^{A,x}
- \tfrac{c_x}{2}(\mu\lambda|\nu\sigma)^{A,x}\Bigr],
$$

$(\cdot)^{A,x}$ は bra が原子 $A$ の殻上で微分されることを意味し、$h^{A,x}$ は勾配の章の 1 次微分
コア Hamiltonian です。動く AO 基底は重なり微分 $S^{A,x}\ne0$ を意味し、これが仮想-占有ブロックが
解かれる*前に*、直交規格性拘束それ自体によって応答の占有-占有ブロックを固定します —
$U^{A,x}_{ij}=-\tfrac12S^{A,x}_{ij}$、$S^{A,x}_{ij}=(C^{\mathsf T}S^{A,x}C)_{ij}$。占有-占有ブロック
の既に既知の応答は、その後、仮想-占有 CPHF 右辺に折り込まれます:

$$
\sum_{bj}(A+B)_{ai,bj}\,U^{A,x}_{bj} = -b^{A,x}_{ai}, \qquad
b^{A,x}_{ai} = F^{A,x}_{ai} - S^{A,x}_{ai}\varepsilon_i - \sum_{kl}(A+B)_{ai,kl}\Bigl(-\tfrac12S^{A,x}_{kl}\Bigr),
$$

これはまさに `linear_solve(apply_aplusb, rhs=-b, ...)` です — 安定性解析と QC-SCF が既に共有する
[同じ CPHF ソルバ](linear-response-cphf.md) であり、この幾何固有の右辺があるだけです。$U^{A,x}$ が
収束すれば、軌道エネルギー応答が副産物として続きます、
$\varepsilon^{A,x}_{ij}=b^{A,x}_{ij}\big|_{oo}+[(A+B)U^{A,x}]_{ij}+U^{A,x}_{ij}(\varepsilon_i-\varepsilon_j)$、
これが折り返しに必要です:

$$
E^{\text{resp}}_{AB} = 4\sum_{\mu\nu}F^{A,x}_{\mu\nu}\bigl(U^{B,y}C_{\text{occ}}^{\mathsf T}\bigr)_{\mu\nu}
- 4\sum_{\mu\nu}W^{B,y}_{\mu\nu}S^{A,x}_{\mu\nu} - 2\sum_{ij}\bigl(C^{\mathsf T}S^{A,x}C\bigr)_{ij}\,
\varepsilon^{B,y}_{ij}.
$$

**ROHF/ROKS は、[勾配の章](analytic-derivatives.md) が $W$ のために導出したのとまったく同じエネルギー
重み付き密度射影子構成を再利用します**(スピン Fock 行列からの $W=\sum_\sigma P_\sigma F_\sigma
P_\sigma$、決して有効 Fock 固有値ではない)— 設計ノートはこれを「[勾配の章] がフラグを立てたのと同じ
微妙さが、今や CPHF 右辺も駆動している」とフラグを立てます、すなわちここで $W$ を間違えると、骨格項
だけでなく結合摂動解法が依存する右辺そのものが壊れます。$W$ の 1 つの誤りが、静かに壊しうる場所は
2 つです。

:::{prf:algorithm} 分子 Hessian の組み立て
:label: alg-hessian-full-ja

**入力:** 収束済み SCF($\mathbf C,\varepsilon,\mathbf D,\mathbf W$)、CPHF 応答エンジン
(`apply_aplusb`)。
**出力:** Hessian $\mathbf H\in\mathbb R^{3n_{\text{atom}}\times3n_{\text{atom}}}$。

1. 骨格ブロック(1 電子+Pulay、3 つの 2 電子 2 次微分クラス、あれば XC)を、固定された
   $\mathbf D$ と $\mathbf W$ と縮約して累積する。
2. すべての原子核 $A$ と Cartesian 方向 $x$ について $F^{A,x}$ と $S^{A,x}$ を構築する(勾配自身の
   `int2e_ip1`/`int1e_ip*` 積分を再利用)。
3. $3n_{\text{atom}}$ 個の摂動それぞれについて:占有-占有ブロック $U^{A,x}_{ij}=-\tfrac12S^{A,x}_{ij}$
   を固定し、それを CPHF 右辺 $b^{A,x}$ に折り込み、共有の `linear_solve` で
   $(A+B)U^{A,x}_{\text{vo}}=-b^{A,x}_{\text{vo}}$ を解く。
4. 収束した $U^{A,x}$ の副産物として軌道エネルギー応答 $\varepsilon^{A,x}$ を計算する。
5. 上記の $E^{\text{resp}}_{AB}$ を通じて各対 $(A,B)$ を Hessian に折り返す。
6. 古典的な核反発 Hessian $\partial^2V_{nn}/\partial R_A\partial R_B$(閉形式)を加える。
7. **自己検定 — 音響和則**:$\sum_B\mathbf H_{AB}\approx\mathbf0$、勾配自身の並進不変性検査の
   Hessian 版(分子全体の剛体並進も 2 次微分を変えない)。
:::

## Hessian から基準振動へ:射影の完全な形

Cartesian Hessian は並進・回転・真の振動を単一の行列に混ぜ込みます — 振動だけを回復するには、単なる
対角化ではなく明示的な射影が必要です。qc-rs の `harmonic_analysis`(`python/qc/thermo.py`、PySCF の
`hessian/thermo.py` の直接移植)は、既製の式ではなく第一原理からこの射影を構成します:

1. **質量重み付け**:$\tilde H_{Ai,Bj}=H_{Ai,Bj}/\sqrt{m_Am_B}$(原子単位の質量)、対称化。
2. **剛体並進/回転基底を質量重み付き座標で直接構築する。** 3 つの並進方向は単に各 Cartesian 軸
   $k$ についての $\sqrt{m_A}\,\hat{\mathbf e}_k$ です;3 つの回転方向は各軸まわりの
   $\sqrt{m_A}\,(\hat{\mathbf e}_k\times\mathbf r_A^{\text{rel}})$ です、$\mathbf r_A^{\text{rel}}$
   は原子の質量中心からの変位です。これら 6 列の QR 分解がそれらを正規直交化し、**線形分子に
   ついては自動的にヌル列を落とします**(分子軸まわりの回転はどの質量も動かさないので、その列は
   厳密にゼロであり QR は消える対角要素を通じてそれを検出します)— これが、qc-rs が $3N-5$ 対
   $3N-6$ を正しくするために「この分子は線形か?」という特殊ケース分岐を決して必要としない*理由*
   です;射影は幾何自体から自動的にそれを扱います。
3. **直交補空間への射影**:$P_{\text{vib}}=I-QQ^{\mathsf T}$、$Q$ は並進/回転基底;$P_{\text{vib}}$
   自体を対角化し固有値 $>0.5$ の固有ベクトルを保持することで、振動部分空間だけの(次元
   $3N-6$ または $3N-5$ の)正規直交基底が得られます。
4. **この振動部分空間に制限された質量重み付き Hessian を対角化する。** 各固有値 $k$ は力の定数です;
   $\tilde\nu=\sqrt{|k|}\cdot(\text{単位因子})$ は波数での振動数を与え、$k$ の*符号*は報告される
   振動数の符号として保持されます — **負の**力の定数は**負の**(すなわち虚数の)報告振動数を与え
   ます、これは真の極小ではなく鞍点の紛れもない特徴です。
5. **質量重みを外す**:生き残った固有ベクトルを実際の Cartesian 原子変位に戻し(再び
   $\sqrt{m_A}$ で割る)、実際に可視化される基準振動を得て、非正規化変位から各モードの**換算質量**
   $\mu_k=1/\sum_{A}|\mathbf d_{A,k}|^2$ を計算します。

## IR 強度:基準振動を直接再利用する

収束済み基準振動と**原子分極率テンソル**(APT、$\partial\mu_i/\partial R_{A,x}$ — 核変位に関する
双極子モーメント微分、RHF/UHF/RKS/UKS と ROHF 参照で利用可能)が与えられれば、モードごとの二重調和
近似 IR 強度は、APT を各質量重み付き基準振動へ直接射影したものです:

$$
\frac{\partial\mu_i}{\partial Q_k} = \sum_{A,x}\frac{\partial\mu_i}{\partial R_{A,x}}\,
\frac{L_{A,x,k}}{\sqrt{m_A}}, \qquad I_k\,[\text{km/mol}] = 974.8801\times4.8032^2\sum_i\Bigl(
\frac{\partial\mu_i}{\partial Q_k}\Bigr)^2,
$$

$L$ は上記射影からの正規直交質量重み付き基準振動行列、定数は原子単位(bohr あたりの電子電荷、
amu$^{-1/2}$ で質量重み付き)を Debye/Å 変換係数経由で慣用的な km/mol 強度単位に変換します。振動数
解析が既に構築したもの以上の新しい機構は不要です — APT はそれ自体、勾配自身の双極子・微分 Fock
積分と密接に関連する微分積分対象であり、強度は APT と基準振動の両方が存在すれば純粋な線形代数縮約
です。

## 理想気体熱化学:厳密な分配関数

調和振動数は標準的な理想気体統計力学に供給されます。qc-rs は CODATA 物理定数(全体を通じて
$k_B$、$h$、$N_A$、Hartree-Joule 変換など)を使い、表を参照するのではなく各分配関数を閉形式で
評価します:

**振動**(調和振動子、実正の各モードにつき 1 因子、振動温度 $\Theta_k=hc\tilde\nu_k/k_B$ と
$x_k=\Theta_k/T$):

$$
\text{ZPE} = \tfrac12\sum_k k_B\Theta_k, \qquad
E_{\text{vib}} = \sum_k k_B\Theta_k\Bigl(\tfrac12+\frac{1}{e^{x_k}-1}\Bigr), \qquad
S_{\text{vib}} = k_B\sum_k\Bigl(\frac{x_k}{e^{x_k}-1}-\ln(1-e^{-x_k})\Bigr).
$$

ここには実正の振動数だけが入ります — 虚数モード(不完全に最適化された幾何から)は無意味な寄与を
生む代わりに振動和から完全に除外されます。これが、`n_imaginary>0` を単に一言添えるべき警告としてで
はなく「これらの熱化学の数値はまだ意味を持たない」として扱うべき理由でもあります。

**並進**(全質量 $M$、温度 $T$、圧力 $P$ の理想気体;これはまさに Sackur-Tetrode の構成です):

$$
q_{\text{trans}} = \Bigl(\frac{2\pi Mk_BT}{h^2}\Bigr)^{3/2}\frac{k_BT}{P}, \qquad
E_{\text{trans}} = \tfrac32k_BT, \qquad S_{\text{trans}} = k_B\bigl(\ln q_{\text{trans}}+\tfrac52\bigr).
$$

**回転**、主慣性モーメント $I_1,I_2,I_3$(同じ核質量/幾何から構築された慣性テンソルの固有値 —
[幾何解析の章](../20-guide/properties/geometric-analysis.md) が回転定数報告のために同じ慣性テンソル
を導出します)と非ゼロモーメントごとの回転温度 $\Theta_{r}=\hbar^2/(2Ik_B)$ から:**線形**分子は
ちょうど 1 つの消える慣性モーメントを持ち(分子軸まわりに回転するエネルギーコストはゼロ)、一般の
非線形式をより単純な式に崩壊させます:

$$
\text{線形:}\quad q_{\text{rot}}=\frac{T}{\sigma\Theta_r}, \quad E_{\text{rot}}=k_BT, \quad
S_{\text{rot}}=k_B\bigl(\ln q_{\text{rot}}+1\bigr);
$$
$$
\text{非線形:}\quad q_{\text{rot}}=\frac{\sqrt\pi}{\sigma}\sqrt{\frac{T^3}{\Theta_{r,1}\Theta_{r,2}
\Theta_{r,3}}}, \quad E_{\text{rot}}=\tfrac32k_BT, \quad S_{\text{rot}}=k_B\bigl(\ln
q_{\text{rot}}+\tfrac32\bigr),
$$

$\sigma$ は回転対称数です(qc-rs は現在 $\sigma=1$ に固定しています — 点群対応の対称数精密化は
既知のフォローアップであり、まだ実装されていません;水のような非対称分子ではこの既定は厳密です
が、対称な分子では補正が必要でしょう)。ある分子が線形と数えられるかどうかは、慣性固有値から
直接決まります(他に比べて消えるほど小さい 1 つ)、上の基準振動射影が使うのと同じ自動検出です —
ユーザー入力から糸を通される別個の「線形か」フラグはありません。

**組み立て**(理想気体関係 $H=U+k_BT$、$G=H-TS$、電子スピン多重度エントロピー項
$S_{\text{elec}}=k_B\ln(\text{mult})$、そして理想気体状態方程式からの $C_p=C_v+k_B$):

$$
U = E_{\text{trans}}+E_{\text{rot}}+E_{\text{vib}}+\text{ZPE}, \qquad H=U+k_BT, \qquad
S = S_{\text{trans}}+S_{\text{rot}}+S_{\text{vib}}+S_{\text{elec}}, \qquad G=H-TS.
$$

検証済み例 — 水/STO-3G、完全な調和+熱化学+IR 強度パイプライン:

```python
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)
np.round(np.sort(fr.frequencies), 2)   # [2041.42, 4494.04, 4796.72]  cm^-1
fr.n_imaginary                          # 0 -- 真の極小
fr.intensities                          # [263.64, 861.86, 459.39]  km/mol(IR 強度)
fr.thermo["zpe"], fr.thermo["g_tot"]     # (0.025817, 0.033348)  Hartree、298.15 K, 101325 Pa
```

3 つの実正振動数、ゼロの虚数モード、そして物理的に妥当な IR 強度の集合が合わせて、これが完全に
自己無撞着な調和解析を伴う真の極小であることを確認します — そして音響和則
($\sum_B H_{AB}\approx0$、[ガイドの章で直接検証](../20-guide/hessian-frequencies-thermo.md))は、
勾配自身の並進不変性検査の Hessian レベルの対応物です。

:::{exercise}
:label: ex-hessian-thermo-theory-ja

1. 並進/回転基底構成は、原子数や幾何への特殊ケース検査ではなく、軸回転列に対して消える対角要素を
   生む QR 分解によって線形分子を検出します。$\hat{\mathbf e}_\parallel$ が線形分子の分子軸のとき、
   $\sqrt{m_A}(\hat{\mathbf e}_\parallel\times\mathbf r_A^{\text{rel}})$ がすべての原子について
   厳密にゼロになる理由、そしてそれがなぜ物理的に線形分子が「回転自由度が 1 つ少ない」ことを意味
   するのかを一文で説明しなさい。
2. 振動エントロピー/エネルギー和は、負(虚数)の振動数を、例えば $|k|$ を取ってともかく含める
   のではなく、明示的に除外します。虚数モードを振動分配関数に含めることが、単に不正確なだけでなく
   物理的に無意味な数値を生む理由は何ですか。
3. 設計ノートは、(勾配を超えて)Hessian に固有の話として、ROHF/ROKS エネルギー重み付き密度
   $W$ を「1 つの誤り、2 つの場所」とフラグ立てます。2 つ目の場所を特定し、そこで $W$ を間違える
   ことが単一のエネルギー項以上を壊す理由を簡潔に説明しなさい。
:::

:::{solution} ex-hessian-thermo-theory-ja
:class: dropdown

1. 線形分子の分子軸まわりに回転しても、すべての原子はちょうど元の位置にとどまります(各原子は
   その軸*上*にあるので、$\mathbf r_A^{\text{rel}}$ は $\hat{\mathbf e}_\parallel$ に平行であり、
   平行な 2 つのベクトルの外積は恒等的にゼロです)— その「回転」に伴う原子変位は、どの原子でも
   どんな質量でも、まったく存在しません。何も変位させない剛体運動は真の自由度ではなく、これが
   まさに「線形分子は回転自由度が 1 つ少ない」ことの意味です:QR 分解の消える対角要素は、単に
   この基底ベクトルが既にゼロベクトルであったことを数値的に検出しているだけです。
2. 調和振動子分配関数とそこから導かれるエントロピー/エネルギー式は、$e^{-\Theta_k/T}$ 型の
   Boltzmann 因子の中に現れる $\Theta_k=hc\tilde\nu_k/k_B$ から構築されます — これらの式は、
   きちんと定義された基底状態と等間隔の励起準位を持つ、真の正曲率調和井戸を仮定しています。負の
   力の定数はそもそも束縛調和振動子に対応しません — それは不安定な方向(鞍点の反応座標)を下る
   運動を記述するものであり、それに対しては分配関数に入れるべききちんと定義された振動基底状態や
   準位間隔が存在しません。$|k|$ を取ってともかく計算すると、エントロピーやエネルギー寄与のように
   見えるが実際にはどの物理的振動モードにも対応しない数値を生みます — これは省略するより悪く、
   幾何が真の極小でないことを(`n_imaginary>0` を通じて)目に見える形で知らせる代わりに、
   $G$ と $S$ を捏造された寄与で静かに汚染してしまうからです。
3. 最初の場所は骨格 Hessian 項自体です(勾配とまったく同様、$W$ は Pulay 重なり微分縮約に入り
   ます)。2 つ目の、Hessian 固有の場所は**CPHF 右辺** $b^{A,x}$ と折り返し項
   $E^{\text{resp}}_{AB}$ であり、どちらも $S^{A,x}$ と縮約された $W^{B,y}$ を明示的に含みます。
   したがって $W$ を間違えることは 1 つの骨格寄与に偏りを与えるだけではありません — $3n_{\text{atom}}$
   回の結合摂動解法すべてを駆動するまさにその摂動を壊すので、誤差は軌道応答 $U^{A,x}$ 自体へ、
   したがってそこから構築されるすべての折り返し項へ伝播します、最終 Hessian の単一の加法的部分
   だけにではありません。
:::

SCF エネルギー、その勾配、その Hessian、そして今や第一原理から導出された完全な振動/熱化学解析を
手にして、[ECP の章](ecp-theory.md) と [溶媒和の章](solvation-theory.md) は、本マニュアルの Part II
がまだ扱っていない残りの Hamiltonian 修正 — 有効内核ポテンシャルと陰的溶媒和 — に向かいます、
どちらも本章が導出した解析的勾配と(実装されている範囲で)解析的 Hessian の機構が既に対応する
ものです。