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

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

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

なぜ 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\) を解く ことで求まります — 線形応答の章 の \(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 グリッドの章で導出 された同じグリッド機構から構築される XC 2 次微分寄与です。

CPHF の右辺と折り返し#

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

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

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

入力: 収束済み 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\)(同じ核質量/幾何から構築された慣性テンソルの固有値 — 幾何解析の章 が回転定数報告のために同じ慣性テンソル を導出します)と非ゼロモーメントごとの回転温度 \(\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 強度パイプライン:

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\)、ガイドの章で直接検証)は、 勾配自身の並進不変性検査の Hessian レベルの対応物です。

練習 12

  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\) を間違える ことが単一のエネルギー項以上を壊す理由を簡潔に説明しなさい。

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