last_modified: 2026-01-08

生成AIによる自動生成記事に関する免責事項: 本記事は、提供された学術論文 The Journal of Chemical Physics, Vol. 49, No. 4, 1719-1729 (1968) および関連する文脈に基づき、大規模言語モデルによって作成された解説記事です。理論的な背景、数学的導出、および数値結果の解釈は、原著論文（J. Gerratt and I. M. Mills）の記述や関連する文献を再構成しています。記事中の数式展開や理論的考察の正確な内容は原著論文を参考にしてください。

1. 序論：量子化学におけるエネルギー微分と応答理論#

分子の分光学的性質や反応性を理解する上で、ポテンシャルエネルギー曲面（PES）の局所的なトポロジー情報は不可欠である。特に、平衡核配置近傍におけるエネルギーの二次微分（ヘシアン、または力の定数）は振動スペクトルの解析に、外部電場に対するエネルギーの一次および二次微分はそれぞれ双極子モーメントおよび分極率に対応する。

1960年代後半、計算機資源の制約の中で、Hartree-Fock (HF) 波動関数を用いてこれらの物性を精度よく算出するための理論的枠組みの構築が急務であった。J. Gerratt と I. M. Mills (1968) は、摂動論的アプローチをHartree-Fock法に適用した Coupled-Perturbed Hartree-Fock (CPHF) 理論を展開し、分子軌道の一次応答（First-order response）を決定する一般化された方程式を導出した 。

本稿では、Gerratt-Millsの定式化に基づき、核変位および外部電場という異なる摂動に対する波動関数の応答を数学的に記述し、それを用いて $H_2$ や LiH などの二原子分子における物性を算出するためのアルゴリズムを詳述する。

---

2. 理論的背景：力の定数とHellmann-Feynman定理の限界#

2.1 分子内の力と力の定数の定義#

$N$ 個の原子核と $n$ 個の電子からなる系を考える。核座標を $R_J (J=1, \dots, N)$、電子座標を $r_\mu (\mu=1, \dots, n)$ とし、Born-Oppenheimer近似の下での電子ハミルトニアン $\mathcal{H}$ を定義する 。 系の全エネルギー $E(R)$ は核座標の関数であり、平衡配置周りでのポテンシャル $V$ は以下のように展開される。

$$V = \frac{1}{2}\sum_{J,K}\sum_{\alpha,\beta} f_{JK}^{\alpha\beta} \delta R_J^\alpha \delta R_K^\beta \tag{1}$$

ここで、$f_{JK}^{\alpha\beta}$ は力の定数（Force Constants）であり、エネルギーの核座標に関する二次微分として定義される 。

$$f_{JK}^{\alpha\beta} = \left( \frac{\partial^2 E}{\partial R_J^\alpha \partial R_K^\beta} \right)_{eq} \tag{2}$$

2.2 ヘルマン-ファインマン項と緩和項#

エネルギーの一次微分については、波動関数 $\Psi$ がハミルトニアンの厳密な固有関数である場合、Hellmann-Feynman定理が成立する。

$$\frac{\partial E}{\partial R_J^\alpha} = \left\langle \Psi \left| \frac{\partial \mathcal{H}}{\partial R_J^\alpha} \right| \Psi \right\rangle \tag{3}$$

これをさらに微分して二次微分（力の定数）を得る際、以下の式が得られる 。

$$f_{JK}^{\alpha\beta} = \frac{\partial^2 \Omega}{\partial R_J^\alpha \partial R_K^\beta} + \int \rho(r, R) \frac{\partial^2 v}{\partial R_J^\alpha \partial R_K^\beta} dr + \int \frac{\partial \rho}{\partial R_K^\beta} \frac{\partial v}{\partial R_J^\alpha} dr \tag{4}$$

ここで、$\Omega$ は核間反発エネルギー、$\rho(r, R)$ は一電子密度関数、$v$ は核ポテンシャルである。 式(4)の右辺第2項までは、電子密度が変化しないと仮定した場合の寄与であり、Hellmann-Feynman項と呼ばれる。一方、第3項は核の動きに伴って電子密度が再分布することによるエネルギーの安定化を表し、**緩和項（Relaxation term）**と呼ばれる 。

対角要素（$J=K$）においては緩和項は常に負の値を取り、化学結合の形成に伴う遮蔽効果などを反映する重要な項である。GerrattとMillsは、Hartree-Fock近似波動関数 $\Phi_{HF}$ を用いた場合でも、この緩和項を正確に見積もるためには、分子軌道（MO）の核座標に対する微分係数 $\partial \phi_i / \partial R$ を知る必要があることを示した 。

---

3. Coupled-Perturbed Hartree-Fock (CPHF) 理論の数理的構造#

CPHF法は、摂動 $\lambda$（核変位や電場）に対するMO係数の応答を決定する手法である。以下にその詳細な導出を行う。

3.1 摂動を含むRoothaan-Hall方程式#

閉殻系におけるHartree-Fock方程式は、基底関数 $\chi_p$ を用いて以下のRoothaan-Hall方程式として記述される 。

$$\mathbf{F}(\lambda) \mathbf{c}_i(\lambda) = \epsilon_i(\lambda) \mathbf{S}(\lambda) \mathbf{c}_i(\lambda) \tag{5}$$ $$\mathbf{c}_i^\dagger(\lambda) \mathbf{S}(\lambda) \mathbf{c}_j(\lambda) = \delta_{ij} \tag{6}$$

ここで、$\lambda$ は摂動パラメータである。$\lambda = \lambda_0$（摂動なし）における解 $\mathbf{c}_i^{(0)}, \epsilon_i^{(0)}$ は既知とする。 問題は、一般の $\lambda$ 近傍におけるMO係数ベクトル $\mathbf{c}_i(\lambda)$ の変化、すなわち一次微分係数を求めることに帰着する 。

3.2 非摂動軌道による展開とユニタリ変換#

計算の簡略化のため、摂動時のMOを、摂動のないMO（Canonical MOs）$\phi_j^{(0)}$ の線形結合として表現する変換を導入する。基底関数での展開係数行列 $\mathbf{C}(\lambda)$ は以下のように書ける 。

$$\mathbf{c}_i(\lambda) = \mathbf{C}^{(0)} \mathbf{u}_i(\lambda) \tag{7}$$

ここで、$\mathbf{u}_i(\lambda)$ は変換後の係数ベクトルである。これにより、方程式は以下のように変換される。

$$\mathfrak{F}(\lambda) \mathbf{u}_i(\lambda) = \epsilon_i(\lambda) \mathbf{O}(\lambda) \mathbf{u}_i(\lambda) \tag{8}$$

ここで、変換されたFock行列 $\mathfrak{F}$ と重なり行列 $\mathbf{O}$ は以下のように定義される。

$$\mathfrak{F}(\lambda) = \mathbf{C}^{(0)\dagger} \mathbf{F}(\lambda) \mathbf{C}^{(0)}, \quad \mathbf{O}(\lambda) = \mathbf{C}^{(0)\dagger} \mathbf{S}(\lambda) \mathbf{C}^{(0)} \tag{9}$$

3.3 摂動展開と一次CPHF方程式#

各項を $\lambda$ についてTaylor展開する 。

$$\mathfrak{F}(\lambda) = \mathfrak{F}^{(0)} + \lambda \mathfrak{F}^{(1)} + \dots \tag{10}$$ $$\mathbf{u}_i(\lambda) = \mathbf{u}_i^{(0)} + \lambda \mathbf{u}_i^{(1)} + \dots \tag{11}$$

$\lambda=0$ において $\mathbf{O}^{(0)} = \mathbf{E}$（単位行列）、$\mathfrak{F}^{(0)}$ は対角行列（対角成分は軌道エネルギー $\epsilon_i^{(0)}$）であることに注意し、一次の項を比較すると以下の関係式が得られる 。

$$(\mathfrak{F}^{(0)} - \epsilon_i^{(0)}\mathbf{E})\mathbf{u}_i^{(1)} = (\epsilon_i^{(0)}\mathbf{O}^{(1)} + \epsilon_i^{(1)}\mathbf{E} - \mathfrak{F}^{(1)})\mathbf{u}_i^{(0)} \tag{12}$$

ここから、混合係数行列の要素 $u_{ij}^{(1)}$ （$j$番目の非摂動MOが$i$番目の摂動MOに混ざる割合）について、以下の重要な式が導かれる 。

$$u_{\sigma j}^{(1)} = \frac{\mathfrak{F}_{\sigma j}^{(1)} - \epsilon_j^{(0)} O_{\sigma j}^{(1)}}{\epsilon_j^{(0)} - \epsilon_\sigma^{(0)}} \quad (\sigma \in \text{virtual}, j \in \text{occupied}) \tag{13}$$

ここで、添字 $\sigma$ は仮想軌道、$j$ は占有軌道を指す。また、対角要素 $u_{jj}^{(1)}$ および占有軌道同士の混合は、正規直交化条件から以下のように決定される 。

$$u_{ij}^{(1)*} + u_{ji}^{(1)} + O_{ji}^{(1)} = 0 \tag{14}$$

3.4 自己無撞着な摂動マトリックスの決定#

式(13)は閉じた式のように見えるが、右辺の $\mathfrak{F}_{\sigma j}^{(1)}$ 自身が一次応答係数 $u_{\tau k}^{(1)}$ に依存しているため、これは連立方程式となる。

Fock行列の一次微分 $\mathfrak{F}_{\sigma j}^{(1)}$ は、一電子部分（コアハミルトニアン）の微分 $\mathcal{H}_{\sigma j}^{(1)}$ と、二電子部分（クーロン・交換積分）の微分 $\mathcal{G}_{\sigma j}^{(1)}$ に分解される 。

$$\mathfrak{F}_{\sigma j}^{(1)} = \mathcal{H}_{\sigma j}^{(1)} + \mathcal{G}_{\sigma j}^{(1)} \tag{15}$$

二電子部分は、密度行列の一次変化（すなわち $u$ 係数）を含むため、以下のような形式になる（詳細は原著論文 Appendix A 参照）。

$$\mathcal{G}_{\sigma j}^{(1)} = \sum_{k=1}^{n} \sum_{\tau=n+1}^{m} \{ 4[\sigma j | k \tau] - [\sigma k | \tau j] - [\sigma \tau | k j] \} u_{\tau k}^{(1)} + \text{(基底関数の微分項)} \tag{16}$$

ここで、記号 $[pq|rs]$ は二電子積分を表す。 最終的に、未知数 $u_{\sigma j}^{(1)}$ に対する線形連立方程式（CPHF方程式）は以下のようにまとめられる 。

$$(\epsilon_j^{(0)} - \epsilon_\sigma^{(0)}) u_{\sigma j}^{(1)} + \sum_{k, \tau} (\text{電子間相互作用項}) u_{\tau k}^{(1)} = \mathcal{H}_{\sigma j}^{(1)} - \epsilon_j^{(0)} O_{\sigma j}^{(1)} + \text{(基底微分項)} \tag{17}$$

この連立方程式を解くことで、波動関数の一次変化が決定される。

---

4. 物理量の計算アルゴリズム：H2およびLiHへの適用#

導出されたCPHF係数を用いて、具体的な物性を計算する手順を示す。ここでは、原著論文で取り上げられている $H_2$ および LiH をモデルケースとする。

4.1 力の定数（Hessian）の計算#

力の定数を求める場合、摂動パラメータ $\lambda$ は核座標 $R_J^\alpha$ である。この場合、基底関数 $\chi_p$ が核中心に固定されているため、基底関数自体も摂動とともに移動する。したがって、積分値の微分には、演算子の微分だけでなく、基底関数の微分 $\partial \chi_p / \partial R$ に由来する項（いわゆるPulay力の一部）が含まれる 。

計算手順:

1. SCF計算: 平衡構造における $\mathbf{C}^{(0)}, \epsilon^{(0)}$ を求める。
2. 積分微分の計算: 一電子積分（$\partial H/\partial R, \partial S/\partial R$）および二電子積分（$\partial [pq|rs]/\partial R$）の一次微分を計算する。これらは解析的に求められる。
3. CPHF方程式の構築: 式(17)の係数行列および右辺ベクトルを構築する。核変位の場合、右辺には積分の微分値が含まれる。
4. 連立方程式の解法: $u_{\sigma j}^{(1)}$ を求める。次元は $(\text{仮想軌道数}) \times (\text{占有軌道数})$ である。
5. エネルギー二次微分の計算: 求まった $u^{(1)}$ を用いて、エネルギーの二次微分式（原著論文 Eq. 54）を評価する 。

$$\frac{\partial E}{\partial \lambda} = 2\sum_{i=1}^n \{ \mathcal{H}_{ii}^{(1)} - \epsilon_i^{(0)}O_{ii}^{(1)} \} + \sum_{i,j=1}^n \frac{\partial}{\partial \lambda} \{ 2[ii|jj] - [ij|ji] \} \tag{18}$$

※ 上記は一次微分（勾配）の式であるが、CPHFの解を用いることで、これをさらに微分した二次微分の表式（緩和項を含む）を評価できる。原著論文では、この緩和項が力の定数の精度に決定的な寄与をすることが示唆されている 。

4.2 双極子モーメント（Dipole Moment）#

双極子モーメント $\mu$ は、ハミルトニアンに対する摂動ではなく、一電子演算子（双極子演算子 $\hat{\mu}$）の期待値として計算される。HF波動関数においては、以下の式で与えられる 。

$$\mu_{HF} = \langle \Phi_{HF} | \hat{\mu} | \Phi_{HF} \rangle = 2 \sum_{i=1}^n \langle \phi_i | \hat{\mu} | \phi_i \rangle \tag{19}$$

これはCPHF方程式を解く必要はなく、収束したSCF波動関数から直接計算される。GerrattとMillsの結果によれば、HFレベルでの双極子モーメントは、真の値に対して一次の誤差（$\epsilon$）ではなく、二次の誤差（$\epsilon^2$）で正しいことが示されており、通常1%程度の精度が得られる 。

4.3 分極率（Polarizability）#

分極率 $\alpha$ は、外部電場 $\mathcal{E}$ に対する双極子モーメントの変化率、すなわちエネルギーの電場に関する二次微分である。

$$\alpha = \frac{\partial \mu}{\partial \mathcal{E}} = \frac{\partial^2 E}{\partial \mathcal{E}^2} \tag{20}$$

計算手順:

1. 摂動パラメータ $\lambda$ を外部電場 $\mathcal{E}$ とする。
2. この場合、基底関数は電場に依存しないため、$\partial \chi_p / \partial \lambda = 0$ となり、重なり行列の微分 $\mathbf{O}^{(1)}$ はゼロになる 。
3. 摂動ハミルトニアンは双極子演算子そのものとなる（$\mathcal{H}^{(1)} = \mathbf{r} \cdot \mathbf{e}$）。
4. 簡略化されたCPHF方程式（原著論文 Eq. 51）を解く 。 $$(\epsilon_j^{(0)} - \epsilon_\sigma^{(0)}) u_{\sigma j}^{(1)} + \sum_{k, \tau} (\dots) u_{\tau k}^{(1)} = \mathcal{H}_{\sigma j}^{(1)} \tag{21}$$
5. 得られた $u^{(1)}$ から分極率を算出する。

4.4 双極子モーメント微分（Dipole Derivatives）#

赤外吸収強度に関係する双極子モーメントの核座標微分 $\partial \mu / \partial R$ は、エネルギーの核座標 $R$ と電場 $\mathcal{E}$ による混合二次微分に対応する。 これは、以下の2通りの方法で計算可能である。

1. 核座標 $R$ に対するCPHFを解き、双極子演算子の期待値の変化を計算する（$\partial \Psi / \partial R$ を使用）。
2. 電場 $\mathcal{E}$ に対するCPHFを解き、Hellmann-Feynman力の電場依存性を計算する（$\partial \Psi / \partial \mathcal{E}$ を使用）。

原著論文では、式(29)として以下の表式を与えている 。

$$\frac{\partial \mu}{\partial R_J^\alpha} = Z_J e_\alpha - 4 \sum_{i=1}^n \left\langle \frac{\partial \phi_i}{\partial R_J^\alpha} \Big| r_\mu \Big| \phi_i \right\rangle \tag{22}$$

ここで、$\partial \phi_i / \partial R$ はCPHFで得られた係数 $u^{(1)}$ と基底関数の微分を含んでいる。

---

5. 結果と考察：Hartree-Fock近似の精度#

GerrattとMillsは、この手法を $H_2$ (水素分子)、LiH (水素化リチウム)、LiF (フッ化リチウム) 等に適用し、実験値との比較を行った。

5.1 $H_2$ における力の定数#

$H_2$ 分子について、異なる基底関数や補間法を用いた結果が比較されている。 表III によれば：

• 実験値: 0.3699 a.u.
• McLean (Curve fitting): 0.4009 a.u.
• Goodisman (Force curve derivative): 0.3799 a.u.

HFレベルでの力の定数は、一般に実験値よりも大きな値（結合が硬い）を与える傾向がある。これは、HF法が電子相関を無視するため、ポテンシャル曲面が実際よりも急峻になる（解離極限での振る舞いが不正確である）ことに起因する。しかし、GerrattとMillsは、HF関数が平衡点近傍では真の波動関数に対する一次の誤差しか持たないため、二次物性（力の定数）も妥当な精度（10%以内）で予測可能であると論じている 。

5.2 緩和項の重要性#

力の定数の計算において、緩和項（Relaxation term）の寄与は無視できない。例えば、LiFのようなイオン性分子においては、核の移動に伴う電子雲の追従（分極）が結合エネルギーに大きく寄与するため、Hellmann-Feynman項のみでは正しくない結果を与えることが示唆されている 。CPHF法は、この電子密度の緩和効果を $u^{(1)}$ を通じて厳密に取り込むことを可能にした。

---

6. 実装に向けたガイドライン#

本稿の読者が実際にプログラムを実装する場合、以下の点に留意する必要がある。

1. 対称性の利用: 分子が対称性を持つ場合、摂動（核変位など）が属する既約表現によって、CPHF方程式はブロック対角化される。これにより計算コストを大幅に削減できる（原著論文 Appendix B 参照） 。
2. 積分微分の実装: Gaussian型軌道（GTO）を用いる現代のプログラムとは異なり、原著ではSlater型軌道（STO）が議論されているが、数理的構造は不変である。積分導関数（Integral Derivatives）の実装が最も煩雑な部分となる。
3. 基底関数の選択: 精度の高い分極率や双極子微分を得るためには、原子価軌道だけでなく、分極関数（Polarization functions）や分散関数（Diffuse functions）を含む柔軟な基底セットが必要である。

---

7. 結論#

GerrattとMillsによる1968年の研究は、Hartree-Fock法に基づき、分子の二次応答物性を解析的に導出するための完全な数理的枠組みを提供した。彼らの導出したCPHF方程式は、現代の量子化学計算ソフトウェア（Gaussian, GAMESS等）において、解析的ヘシアン（Analytic Hessian）やNMR遮蔽定数などの計算エンジンの基礎となっている。

この理論は、単に数値を出すだけでなく、「核が動いたときに電子雲がどのように歪むか（Relaxation）」、「外部電場に対して軌道がどう混合するか（Polarization）」といった物理的直観を、厳密な数学的言葉で記述することを可能にした点で、理論化学における重要なマイルストーンである。

参考文献#

本記事は、以下の文献の記述に基づき作成された。

• [1] J. Gerratt and I. M. Mills, “Force Constants and Dipole-Moment Derivatives of Molecules from Perturbed Hartree-Fock Calculations. I”, J. Chem. Phys. 49, 1719-1729 (1968).

• [2] R. McWeeny, Methods of Molecular Quantum Mechanics, Academic Press, London (1992).

• [3] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory, Dover Publications, New York (1996).

• [4] C. A. Coulson and H. C. Longuet-Higgins, “The electronic structure of conjugated systems. I. General theory”, Proc. Roy. Soc. A 191, 39–60 (1947).

• [5] P. Pulay, “Ab initio calculation of force constants and equilibrium geometries in polyatomic molecules”, Mol. Phys. 17, 197–204 (1969).

• [6] P. Pulay, “Analytical derivative methods in quantum chemistry”, Adv. Chem. Phys. 69, 241–293 (1987).

• [7] N. C. Handy and H. F. Schaefer, “On the evaluation of analytic energy derivatives for correlated wave functions”, J. Chem. Phys. 81, 5031–5033 (1984).

• [8] D. M. Bishop, Group Theory and Chemistry, Dover Publications, New York (1993).

• [9] C. E. Dykstra, G. Frenking, K. S. Kim, and G. E. Scuseria (eds.), Theory and Applications of Computational Chemistry, Elsevier, Amsterdam (2005).

• [10] M. Karplus and H. J. Kolker, “Response of molecular systems to external perturbations”, J. Chem. Phys. 41, 3955–3970 (1964).

• [11] J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Binkley, “Derivative studies in Hartree–Fock and Møller–Plesset theories”, Int. J. Quantum Chem. 13, 225–241 (1978).

• [12] Y. Yamaguchi, Y. Osamura, J. D. Goddard III, and H. F. Schaefer III, A New Dimension to Quantum Chemistry: Analytic Derivative Methods, Oxford University Press, Oxford (1994).

• [13] T. Helgaker, P. Jørgensen, and J. Olsen, Molecular Electronic-Structure Theory, Wiley, Chichester (2000).

• [14] M. Head-Gordon and J. A. Pople, “A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations”, J. Chem. Phys. 89, 5777–5786 (1988).

• [15] C. Ochsenfeld, J. Kussmann, and D. S. Lambrecht, “Linear-scaling techniques in ab initio quantum chemistry”, Chem. Rev. 107, 3589–3614 (2007).

補足#

cphf.py

import numpy as np
from scipy.special import erf
from scipy.linalg import solve, LinAlgError

# =============================================================================
# 1. Mathematical Utilities
# =============================================================================

def F_n(n, t):
    """
    Boys Function F_n(t) = integral_0^1 u^(2n) exp(-t u^2) du.
    Using asymptotic expansion for small t to avoid numerical instability.
    """
    if t < 1e-8:
        return 1.0 / (2*n + 1) - t / (2*n + 3)
    else:
        f0 = 0.5 * (np.pi/t)**0.5 * erf(t**0.5)
        if n == 0: return f0
        
        # F_{n+1} = [(2n+1)F_n - exp(-t)] / (2t)
        exp_t = np.exp(-t)
        f1 = (f0 - exp_t) / (2*t)
        if n == 1: return f1
        
        f2 = (3*f1 - exp_t) / (2*t)
        if n == 2: return f2
        
        raise ValueError("Order n > 2 not implemented.")

def get_F0_F1_F2(t):
    """Helper to return F0, F1, F2 simultaneously."""
    if t < 1e-8:
        f0 = 1.0 - t/3.0
        f1 = 1.0/3.0 - t/5.0
        f2 = 1.0/5.0 - t/7.0
        return f0, f1, f2
    
    val_sqrt = t**0.5
    exp_t = np.exp(-t)
    f0 = 0.5 * (np.pi/t)**0.5 * erf(val_sqrt)
    f1 = (f0 - exp_t) / (2*t)
    f2 = (3*f1 - exp_t) / (2*t)
    return f0, f1, f2

# =============================================================================
# 2. Molecular Objects
# =============================================================================

class BasisFunction:
    def __init__(self, origin, params):
        self.origin = np.array(origin, dtype=float)
        self.exps = np.array(params[0])
        self.coefs = np.array(params[1])
        
        # Pre-compute normalized coefficients
        norm_coefs_list = []
        for a, c in zip(self.exps, self.coefs):
            N = (2*a/np.pi)**0.75
            norm_coefs_list.append(c * N)
        self.norm_coefs = np.array(norm_coefs_list) # Convert to numpy array for speed

class Atom:
    def __init__(self, symbol, coords, atomic_num):
        self.symbol = symbol
        self.coords = np.array(coords, dtype=float)
        self.Z = atomic_num
        self.basis = []

# =============================================================================
# 3. Integral Engine (0th Derivative)
# =============================================================================

class IntegralEngine:
    def __init__(self, atoms, basis_flat):
        self.atoms = atoms
        self.basis = basis_flat
        self.n_basis = len(basis_flat)

    def compute_integrals(self):
        S = np.zeros((self.n_basis, self.n_basis))
        T = np.zeros((self.n_basis, self.n_basis))
        V = np.zeros((self.n_basis, self.n_basis))
        ERI = np.zeros((self.n_basis, self.n_basis, self.n_basis, self.n_basis))

        for i in range(self.n_basis):
            for j in range(i, self.n_basis):
                s_val, t_val, v_val = self._overlap_kinetic_nuc(self.basis[i], self.basis[j])
                S[i, j] = S[j, i] = s_val
                T[i, j] = T[j, i] = t_val
                V[i, j] = V[j, i] = v_val

        # Naive O(N^4) loop (Symmetry could reduce this by ~8x)
        for i in range(self.n_basis):
            for j in range(self.n_basis):
                for k in range(self.n_basis):
                    for l in range(self.n_basis):
                        ERI[i, j, k, l] = self._eri_quartet(self.basis[i], self.basis[j], 
                                                            self.basis[k], self.basis[l])
        return S, T, V, ERI

    def _overlap_kinetic_nuc(self, bi, bj):
        s_val = t_val = v_val = 0.0
        vec_AB = bi.origin - bj.origin
        dist_sq = np.dot(vec_AB, vec_AB)

        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b
                mu = a * b / p
                rp = (a*bi.origin + b*bj.origin) / p
                K_ab = np.exp(-mu * dist_sq)
                
                ov = (np.pi/p)**1.5 * K_ab
                kin = mu * (3.0 - 2.0*mu*dist_sq) * ov
                
                nuc = 0.0
                for atom in self.atoms:
                    vec_PC = rp - atom.coords
                    T_val = p * np.dot(vec_PC, vec_PC)
                    nuc += -atom.Z * (2*np.pi/p) * K_ab * F_n(0, T_val)
                
                term = ca * cb
                s_val += term * ov
                t_val += term * kin
                v_val += term * nuc
        return s_val, t_val, v_val

    def _eri_quartet(self, bi, bj, bk, bl):
        val = 0.0
        vec_AB = bi.origin - bj.origin
        vec_CD = bk.origin - bl.origin
        R_AB2 = np.dot(vec_AB, vec_AB)
        R_CD2 = np.dot(vec_CD, vec_CD)

        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                p = bi.exps[ia] + bj.exps[jb]
                rp = (bi.exps[ia]*bi.origin + bj.exps[jb]*bj.origin) / p
                # Precompute K_ab part
                K_ab = np.exp(-(bi.exps[ia]*bj.exps[jb]/p) * R_AB2)
                
                for kc, cc in enumerate(bk.norm_coefs):
                    for ld, cd in enumerate(bl.norm_coefs):
                        q = bk.exps[kc] + bl.exps[ld]
                        rq = (bk.exps[kc]*bk.origin + bl.exps[ld]*bl.origin) / q
                        K_cd = np.exp(-(bk.exps[kc]*bl.exps[ld]/q) * R_CD2)

                        dist_PQ2 = np.linalg.norm(rp - rq)**2
                        alpha = p * q / (p + q)
                        T_val = alpha * dist_PQ2
                        
                        term = (2 * np.pi**2.5 / (p*q*(p+q)**0.5)) * F_n(0, T_val)
                        val += ca * cb * cc * cd * K_ab * K_cd * term
        return val

# =============================================================================
# 4. Analytic 1st Derivatives (Gradient)
# =============================================================================

class AnalyticGradient:
    def __init__(self, atoms, basis_flat):
        self.atoms = atoms
        self.basis = basis_flat
        self.n_basis = len(basis_flat)
        self.n_atoms = len(atoms)
        
        # Strict mapping check
        self.basis_atom_map = []
        for bf in basis_flat:
            mapped = False
            for i, at in enumerate(atoms):
                if np.linalg.norm(bf.origin - at.coords) < 1e-6:
                    self.basis_atom_map.append(i)
                    mapped = True
                    break
            if not mapped:
                raise ValueError(f"Basis function at {bf.origin} not mapped to any atom.")

    def run(self):
        n_coords = 3 * self.n_atoms
        dS = np.zeros((n_coords, self.n_basis, self.n_basis))
        dH = np.zeros((n_coords, self.n_basis, self.n_basis))
        dERI = np.zeros((n_coords, self.n_basis, self.n_basis, self.n_basis, self.n_basis))

        # --- 1-Electron Integrals (Fix: Prevent double counting for diagonal) ---
        for i in range(self.n_basis):
            for j in range(i, self.n_basis): # Loop upper triangle
                bi, bj = self.basis[i], self.basis[j]
                idx_A, idx_B = self.basis_atom_map[i], self.basis_atom_map[j]
                
                # Get gradients for S, T, V
                ds_dict, dt_dict, dv_dict = self._grad_1e(bi, bj, idx_A, idx_B)
                
                # Merge T and V into H
                dh_dict = {}
                for k in dt_dict: dh_dict[k] = dt_dict[k]
                for k in dv_dict: 
                    if k in dh_dict: dh_dict[k] += dv_dict[k]
                    else: dh_dict[k] = dv_dict[k]
                
                # Fill Tensors
                for atom_idx, val_vec in ds_dict.items():
                    for dim in range(3):
                        row = atom_idx * 3 + dim
                        val = val_vec[dim]
                        dS[row, i, j] += val
                        if i != j:         
                            dS[row, j, i] += val
                            
                for atom_idx, val_vec in dh_dict.items():
                    for dim in range(3):
                        row = atom_idx * 3 + dim
                        val = val_vec[dim]
                        dH[row, i, j] += val
                        if i != j:         
                            dH[row, j, i] += val

        # --- 2-Electron Integrals ---
        # (Assuming original implementation loops over necessary quartets.
        #  If original uses i,j,k,l full loops, it's fine.
        #  If it uses upper triangle loops like i<=j, k<=l, similar checks are needed.)
        #  For safety in this context, assuming standard full loop or managed quartets:
        for i in range(self.n_basis):
            for j in range(self.n_basis):
                for k in range(self.n_basis):
                    for l in range(self.n_basis):
                        grads = self._grad_eri_quartet(self.basis[i], self.basis[j], 
                                                       self.basis[k], self.basis[l])
                        for atom_idx, g_vec in grads.items():
                            for ax in range(3):
                                dERI[atom_idx*3+ax, i, j, k, l] += g_vec[ax]
                                
        return dS, dH, dERI

    def _grad_1e(self, bi, bj, idx_A, idx_B):
        ds = {idx_A: np.zeros(3), idx_B: np.zeros(3)}
        dt = {idx_A: np.zeros(3), idx_B: np.zeros(3)}
        dv = {} 
        
        A, B = bi.origin, bj.origin
        vec_AB = A - B
        R2 = np.dot(vec_AB, vec_AB)

        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b
                xi = a * b / p
                rp = (a*A + b*B) / p
                K = np.exp(-xi * R2)
                
                S_prim = (np.pi/p)**1.5 * K
                dS_dA = -2 * xi * vec_AB * S_prim
                
                term1 = xi * (3 - 2*xi*R2)
                d_term1 = -4 * xi**2 * vec_AB
                dT_dA = d_term1 * S_prim + term1 * dS_dA
                
                term = ca * cb
                ds[idx_A] += term * dS_dA
                ds[idx_B] -= term * dS_dA
                
                dt[idx_A] += term * dT_dA
                dt[idx_B] -= term * dT_dA
                
                for k_at, atom in enumerate(self.atoms):
                    C = atom.coords
                    Z = atom.Z
                    vec_PC = rp - C
                    T = p * np.dot(vec_PC, vec_PC)
                    Pre = -Z * (2*np.pi/p)
                    f0, f1, _ = get_F0_F1_F2(T)
                    
                    dK_dA = -2 * xi * vec_AB * K
                    dTheta_dA = 2 * a * vec_PC
                    dV_dA = Pre * (dK_dA * f0 - K * f1 * dTheta_dA)
                    
                    dK_dB = 2 * xi * vec_AB * K
                    dTheta_dB = 2 * b * vec_PC
                    dV_dB = Pre * (dK_dB * f0 - K * f1 * dTheta_dB)
                    
                    dTheta_dC = -2 * p * vec_PC
                    dV_dC = Pre * K * (-f1) * dTheta_dC
                    
                    if idx_A not in dv: dv[idx_A] = np.zeros(3)
                    dv[idx_A] += term * dV_dA
                    if idx_B not in dv: dv[idx_B] = np.zeros(3)
                    dv[idx_B] += term * dV_dB
                    if k_at not in dv: dv[k_at] = np.zeros(3)
                    dv[k_at] += term * dV_dC
                    
        return ds, dt, dv

    def _grad_eri_quartet(self, bi, bj, bk, bl):
        idx_map = [self.basis_atom_map[self.basis.index(x)] for x in [bi, bj, bk, bl]]
        grads = {i: np.zeros(3) for i in idx_map}
        
        A, B = bi.origin, bj.origin
        C, D = bk.origin, bl.origin
        vec_AB, vec_CD = A - B, C - D
        R_AB2, R_CD2 = np.dot(vec_AB, vec_AB), np.dot(vec_CD, vec_CD)

        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                p = bi.exps[ia] + bj.exps[jb]
                rp = (bi.exps[ia]*A + bj.exps[jb]*B) / p
                xi_ab = bi.exps[ia]*bj.exps[jb]/p
                K_ab = np.exp(-xi_ab * R_AB2)
                dK_dA = -2 * xi_ab * vec_AB * K_ab

                for kc, cc in enumerate(bk.norm_coefs):
                    for ld, cd in enumerate(bl.norm_coefs):
                        q = bk.exps[kc] + bl.exps[ld]
                        rq = (bk.exps[kc]*C + bl.exps[ld]*D) / q
                        xi_cd = bk.exps[kc]*bl.exps[ld]/q
                        K_cd = np.exp(-xi_cd * R_CD2)
                        dK_dC = -2 * xi_cd * vec_CD * K_cd
                        
                        vec_PQ = rp - rq
                        alpha = p * q / (p + q)
                        T = alpha * np.dot(vec_PQ, vec_PQ)
                        
                        f0, f1, _ = get_F0_F1_F2(T)
                        coef = ca * cb * cc * cd
                        Pre = 2 * np.pi**2.5 / (p*q*(p+q)**0.5)
                        
                        dT_dP = 2 * alpha * vec_PQ
                        dT_dQ = -2 * alpha * vec_PQ
                        
                        term_A = dK_dA * K_cd * f0 - K_ab * K_cd * f1 * dT_dP * (bi.exps[ia]/p)
                        grads[idx_map[0]] += coef * Pre * term_A
                        
                        term_B = -dK_dA * K_cd * f0 - K_ab * K_cd * f1 * dT_dP * (bj.exps[jb]/p)
                        grads[idx_map[1]] += coef * Pre * term_B
                        
                        term_C = K_ab * dK_dC * f0 - K_ab * K_cd * f1 * dT_dQ * (bk.exps[kc]/q)
                        grads[idx_map[2]] += coef * Pre * term_C
                        
                        term_D = K_ab * (-dK_dC) * f0 - K_ab * K_cd * f1 * dT_dQ * (bl.exps[ld]/q)
                        grads[idx_map[3]] += coef * Pre * term_D
        return grads

# =============================================================================
# 5. Analytic Hessian (2nd Derivative)
# =============================================================================

class AnalyticHessian:
    def __init__(self, atoms, basis_flat):
        self.atoms = atoms
        self.basis = basis_flat
        self.n_basis = len(basis_flat)
        self.n_atoms = len(atoms)
        self.n_coords = 3 * self.n_atoms
        self.basis_atom_map = []
        for bf in basis_flat:
            mapped = False
            for i, at in enumerate(atoms):
                if np.linalg.norm(bf.origin - at.coords) < 1e-6:
                    self.basis_atom_map.append(i)
                    mapped = True
                    break
            if not mapped:
                raise ValueError(f"Basis function at {bf.origin} not mapped to any atom.")
    
    def run(self):
        print("Computing Analytic Hessian...")
        d2S = np.zeros((self.n_coords, self.n_coords, self.n_basis, self.n_basis))
        d2H = np.zeros((self.n_coords, self.n_coords, self.n_basis, self.n_basis))
        
        # 1-Electron
        for i in range(self.n_basis):
            for j in range(i, self.n_basis):
                bi, bj = self.basis[i], self.basis[j]
                idx_A, idx_B = self.basis_atom_map[i], self.basis_atom_map[j]
                
                s_dict = self._hess_S(bi, bj, idx_A, idx_B)
                h_dict = self._hess_T(bi, bj, idx_A, idx_B)
                v_dict = self._hess_V(bi, bj, idx_A, idx_B)
                
                # Merge V into H
                for key, val in v_dict.items():
                    if key in h_dict: h_dict[key] += val
                    else: h_dict[key] = val
                
                # Fill Tensors (Symmetric ij)
                for (u, v), mat in s_dict.items():
                    r1, r2 = u*3, v*3
                    for a in range(3):
                        for b in range(3):
                            d2S[r1+a, r2+b, i, j] += mat[a,b]
                            d2S[r1+a, r2+b, j, i] += mat[a,b]
                            
                for (u, v), mat in h_dict.items():
                    r1, r2 = u*3, v*3
                    for a in range(3):
                        for b in range(3):
                            d2H[r1+a, r2+b, i, j] += mat[a,b]
                            d2H[r1+a, r2+b, j, i] += mat[a,b]

        # 2-Electron
        d2ERI = self._hess_eri()
        return d2S, d2H, d2ERI

    def _hess_S(self, bi, bj, idx_A, idx_B):
        res = {}
        A, B = bi.origin, bj.origin
        vec_AB = A - B
        R2 = np.dot(vec_AB, vec_AB)
        I3 = np.eye(3)
        Q_outer = np.outer(vec_AB, vec_AB)
        mat_AA = np.zeros((3,3))
        
        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b
                xi = a * b / p
                S_prim = (np.pi/p)**1.5 * np.exp(-xi * R2)
                mat_AA += ca * cb * S_prim * (4 * xi**2 * Q_outer - 2 * xi * I3)
        
        res[(idx_A, idx_A)] = mat_AA
        res[(idx_B, idx_B)] = mat_AA
        res[(idx_A, idx_B)] = -mat_AA
        res[(idx_B, idx_A)] = -mat_AA
        return res

    def _hess_T(self, bi, bj, idx_A, idx_B):
        res = {}
        A, B = bi.origin, bj.origin
        vec_AB = A - B
        R2 = np.dot(vec_AB, vec_AB)
        I3 = np.eye(3)
        Q_outer = np.outer(vec_AB, vec_AB)
        mat_AA = np.zeros((3,3))
        
        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b
                xi = a * b / p
                S_prim = (np.pi/p)**1.5 * np.exp(-xi * R2)
                
                F = 3*xi - 2*xi**2 * R2
                dF_dQ = -4 * xi**2 * vec_AB
                d2F_dQ2 = -4 * xi**2 * I3
                dS_dQ = -2 * xi * vec_AB * S_prim
                d2S_dQ2 = S_prim * (4 * xi**2 * Q_outer - 2 * xi * I3)
                
                term = (d2F_dQ2 * S_prim) + \
                       (np.outer(dF_dQ, dS_dQ) + np.outer(dS_dQ, dF_dQ)) + \
                       (F * d2S_dQ2)
                mat_AA += ca * cb * term
                
        res[(idx_A, idx_A)] = mat_AA
        res[(idx_B, idx_B)] = mat_AA
        res[(idx_A, idx_B)] = -mat_AA
        res[(idx_B, idx_A)] = -mat_AA
        return res

    def _hess_V(self, bi, bj, idx_A, idx_B):
        res = {}
        A, B = bi.origin, bj.origin
        vec_AB = A - B
        R2 = np.dot(vec_AB, vec_AB)
        I3 = np.eye(3)
        
        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b
                rp = (a*A + b*B) / p
                xi = a * b / p
                K = np.exp(-xi * R2)
                
                dK_dA = -2 * xi * vec_AB * K
                d2K_AA = K * (4 * xi**2 * np.outer(vec_AB, vec_AB) - 2 * xi * I3)
                
                for k_at, atom in enumerate(self.atoms):
                    C = atom.coords
                    Z = atom.Z
                    idx_C = k_at
                    
                    vec_PC = rp - C
                    T_val = p * np.dot(vec_PC, vec_PC)
                    Pre = -Z * (2*np.pi/p)
                    f0, f1, f2 = get_F0_F1_F2(T_val)
                    
                    wa, wb = a/p, b/p
                    dT_dA = 2 * p * wa * vec_PC
                    dT_dB = 2 * p * wb * vec_PC
                    dT_dC = -2 * p * vec_PC
                    
                    d2T_AA = 2 * p * wa * wa * I3
                    d2T_BB = 2 * p * wb * wb * I3
                    d2T_CC = 2 * p * I3
                    d2T_AB = 2 * p * wa * wb * I3
                    d2T_AC = -2 * p * wa * I3
                    d2T_BC = -2 * p * wb * I3
                    
                    indices = [idx_A, idx_B, idx_C]
                    unique = sorted(list(set(indices)))
                    
                    for u in unique:
                        for v in unique:
                            k_uv = np.zeros((3,3))
                            if (u!=idx_C) and (v!=idx_C):
                                k_uv = d2K_AA if u==v else -d2K_AA
                            
                            k_u = np.zeros(3)
                            if u==idx_A: k_u = dK_dA
                            elif u==idx_B: k_u = -dK_dA
                            
                            k_v = np.zeros(3)
                            if v==idx_A: k_v = dK_dA
                            elif v==idx_B: k_v = -dK_dA
                            
                            t_u = dT_dA if u==idx_A else (dT_dB if u==idx_B else dT_dC)
                            t_v = dT_dA if v==idx_A else (dT_dB if v==idx_B else dT_dC)
                            
                            t_uv = np.zeros((3,3))
                            if u==idx_A and v==idx_A: t_uv = d2T_AA
                            elif u==idx_B and v==idx_B: t_uv = d2T_BB
                            elif u==idx_C and v==idx_C: t_uv = d2T_CC
                            elif (u,v) in [(idx_A,idx_B),(idx_B,idx_A)]: t_uv = d2T_AB
                            elif (u,v) in [(idx_A,idx_C),(idx_C,idx_A)]: t_uv = d2T_AC
                            elif (u,v) in [(idx_B,idx_C),(idx_C,idx_B)]: t_uv = d2T_BC
                            
                            term = k_uv * f0
                            term += np.outer(k_u, t_v) * (-f1)
                            term += np.outer(t_u, k_v) * (-f1)
                            term += K * (f2 * np.outer(t_u, t_v) - f1 * t_uv)
                            
                            if (u,v) not in res: res[(u,v)] = np.zeros((3,3))
                            res[(u,v)] += ca * cb * Pre * term
        return res

    def _hess_eri(self):
        d2ERI = np.zeros((self.n_coords, self.n_coords, self.n_basis, self.n_basis, self.n_basis, self.n_basis))
        n = self.n_basis
        for i in range(n):
            for j in range(n):
                for k in range(n):
                    for l in range(n):
                        grads = self._hess_eri_quartet(self.basis[i], self.basis[j], 
                                                       self.basis[k], self.basis[l])
                        for (u, v), mat in grads.items():
                            r1, r2 = u*3, v*3
                            for a in range(3):
                                for b in range(3):
                                    d2ERI[r1+a, r2+b, i, j, k, l] += mat[a,b]
        return d2ERI

    def _hess_eri_quartet(self, bi, bj, bk, bl):
        res = {}
        idx_map = [self.basis_atom_map[self.basis.index(x)] for x in [bi, bj, bk, bl]]
        idx_A, idx_B, idx_C, idx_D = idx_map
        unique = sorted(list(set(idx_map)))
        
        A, B = bi.origin, bj.origin
        C, D = bk.origin, bl.origin
        vec_AB, vec_CD = A-B, C-D
        R_AB2, R_CD2 = np.dot(vec_AB, vec_AB), np.dot(vec_CD, vec_CD)
        I3 = np.eye(3)

        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                p = bi.exps[ia] + bj.exps[jb]
                xi_ab = bi.exps[ia]*bj.exps[jb]/p
                K_ab = np.exp(-xi_ab * R_AB2)
                dk1_A = -2*xi_ab*vec_AB*K_ab
                d2k1_AA = K_ab * (4*xi_ab**2 * np.outer(vec_AB, vec_AB) - 2*xi_ab*I3)
                
                for kc, cc in enumerate(bk.norm_coefs):
                    for ld, cd in enumerate(bl.norm_coefs):
                        q = bk.exps[kc] + bl.exps[ld]
                        xi_cd = bk.exps[kc]*bl.exps[ld]/q
                        K_cd = np.exp(-xi_cd * R_CD2)
                        dk2_C = -2*xi_cd*vec_CD*K_cd
                        d2k2_CC = K_cd * (4*xi_cd**2 * np.outer(vec_CD, vec_CD) - 2*xi_cd*I3)
                        
                        rp = (bi.exps[ia]*A + bj.exps[jb]*B)/p
                        rq = (bk.exps[kc]*C + bl.exps[ld]*D)/q
                        vec_PQ = rp - rq
                        alpha = p*q/(p+q)
                        T_val = alpha * np.dot(vec_PQ, vec_PQ)
                        
                        f0, f1, f2 = get_F0_F1_F2(T_val)
                        coef = ca * cb * cc * cd
                        Pre = 2 * np.pi**2.5 / (p*q*(p+q)**0.5)
                        
                        wa, wb = bi.exps[ia]/p, bj.exps[jb]/p
                        wc, wd = bk.exps[kc]/q, bl.exps[ld]/q
                        
                        for u in unique:
                            for v in unique:
                                k1_u = dk1_A if u==idx_A else (-dk1_A if u==idx_B else np.zeros(3))
                                k1_v = dk1_A if v==idx_A else (-dk1_A if v==idx_B else np.zeros(3))
                                k2_u = dk2_C if u==idx_C else (-dk2_C if u==idx_D else np.zeros(3))
                                k2_v = dk2_C if v==idx_C else (-dk2_C if v==idx_D else np.zeros(3))
                                
                                k1_uv = np.zeros((3,3))
                                if (u in [idx_A, idx_B]) and (v in [idx_A, idx_B]):
                                    sign = 1 if u==v else -1
                                    k1_uv = sign * d2k1_AA
                                    
                                k2_uv = np.zeros((3,3))
                                if (u in [idx_C, idx_D]) and (v in [idx_C, idx_D]):
                                    sign = 1 if u==v else -1
                                    k2_uv = sign * d2k2_CC

                                def get_w(idx):
                                    w = 0
                                    if idx==idx_A: w += wa
                                    if idx==idx_B: w += wb
                                    if idx==idx_C: w -= wc
                                    if idx==idx_D: w -= wd
                                    return w
                                
                                w_u, w_v = get_w(u), get_w(v)
                                t_u = 2 * alpha * vec_PQ * w_u
                                t_v = 2 * alpha * vec_PQ * w_v
                                t_uv = 2 * alpha * w_u * w_v * I3
                                
                                term = (k1_uv * K_cd + K_ab * k2_uv + np.outer(k1_u, k2_v) + np.outer(k2_u, k1_v)) * f0
                                term += np.outer(k1_u * K_cd + K_ab * k2_u, t_v) * (-f1)
                                term += np.outer(t_u, k1_v * K_cd + K_ab * k2_v) * (-f1)
                                term += K_ab * K_cd * (f2 * np.outer(t_u, t_v) - f1 * t_uv)
                                
                                if (u,v) not in res: res[(u,v)] = np.zeros((3,3))
                                res[(u,v)] += coef * Pre * term
        return res

# =============================================================================
# 6. CPHF Solver & Hessian Driver
# =============================================================================

def run_calculation(molecule_name, r_bond, verify=False):
    print(f"\n=== Calculation for {molecule_name} (R={r_bond}) ===")
    
    # 1. Setup Molecule
    if molecule_name == 'H2':
        h1 = Atom('H', [0,0,0], 1)
        h2 = Atom('H', [0,0,r_bond], 1)
        exps = [3.42525, 0.623913, 0.168855]
        coefs = [0.154329, 0.535328, 0.444635]
        h1.basis.append(BasisFunction(h1.coords, (exps, coefs)))
        h2.basis.append(BasisFunction(h2.coords, (exps, coefs)))
        atoms = [h1, h2]
    elif molecule_name == 'LiH':
        li = Atom('Li', [0,0,0], 3)
        h = Atom('H', [0,0,r_bond], 1)
        # Using H-like 1s for simplification in demo
        exps_li = [16.119575, 2.9362007, 0.7946505]
        coefs_li = [0.15432897, 0.53532814, 0.44463454]
        li.basis.append(BasisFunction(li.coords, (exps_li, coefs_li)))
        
        exps_h = [3.42525, 0.623913, 0.168855]
        coefs_h = [0.154329, 0.535328, 0.444635]
        h.basis.append(BasisFunction(h.coords, (exps_h, coefs_h)))
        atoms = [li, h]
    
    basis_flat = [b for a in atoms for b in a.basis]
    n_basis = len(basis_flat)
    n_occ = sum([a.Z for a in atoms]) // 2
    n_vir = n_basis - n_occ

    # 2. SCF (RHF)
    eng = IntegralEngine(atoms, basis_flat)
    S, T, V, ERI = eng.compute_integrals()
    Hcore = T + V
    
    # --- FIX: Numerical stability for Orthogonalizer ---
    evals, evecs = np.linalg.eigh(S)
    eps_thresh = 1e-8
    evals = np.clip(evals, eps_thresh, None) # Clip small eigenvalues
    X = evecs @ np.diag(1.0/np.sqrt(evals)) @ evecs.T
    
    P = np.zeros((n_basis, n_basis))
    C = np.zeros((n_basis, n_basis))
    eps = np.zeros(n_basis)
    E_elec = 0.0
    
    for itr in range(50):
        G = np.einsum('uvls,ls->uv', ERI, P)
        K_mat = np.einsum('ulvs,ls->uv', ERI, P)
        F = Hcore + 2*G - K_mat
        
        F_p = X.T @ F @ X
        eps, C_p = np.linalg.eigh(F_p)
        C = X @ C_p
        C_occ = C[:, :n_occ]
        P_new = C_occ @ C_occ.T
        
        E_new = np.sum(P_new * (Hcore + F)) * 0.5
        if abs(E_new - E_elec) < 1e-9: break
        P = P_new
        E_elec = E_new
    
    E_nuc = 0.0
    for i in range(len(atoms)):
        for j in range(i+1, len(atoms)):
            dist = np.linalg.norm(atoms[i].coords - atoms[j].coords)
            E_nuc += atoms[i].Z * atoms[j].Z / dist
    
    print(f"  E_HF = {E_elec + E_nuc:.6f} Hartree")

    # 3. Analytic Gradient
    ag = AnalyticGradient(atoms, basis_flat)
    dS, dH, dERI = ag.run()

    # 4. Analytic Hessian
    ah = AnalyticHessian(atoms, basis_flat)
    d2S, d2H, d2ERI = ah.run()

    # 5. CPHF Solver
    tmp = np.einsum('uvls, up->pvls', ERI, C)
    tmp = np.einsum('pvls, vq->pqls', tmp, C)
    tmp = np.einsum('pqls, lr->pqrs', tmp, C)
    I_mo = np.einsum('pqrs, st->pqrt', tmp, C)
    
    dim = n_occ * n_vir
    A_mat = np.zeros((dim, dim))
    map_ai = []
    for a in range(n_vir):
        for i in range(n_occ):
            map_ai.append((a, i))
            
    for idx1 in range(dim):
        a, i = map_ai[idx1]
        a_mo, i_mo = n_occ + a, i
        for idx2 in range(dim):
            b, j = map_ai[idx2]
            b_mo, j_mo = n_occ + b, j
            val = 0.0
            if idx1 == idx2: val += eps[a_mo] - eps[i_mo]
            val += 4 * I_mo[a_mo, i_mo, b_mo, j_mo]
            val -= I_mo[a_mo, b_mo, j_mo, i_mo]
            val -= I_mo[a_mo, j_mo, b_mo, i_mo]
            A_mat[idx1, idx2] = val
    
    # --- FIX: Use solve instead of inv + regularization fallback ---
    U = np.zeros((3*len(atoms), n_vir, n_occ))
    
    for lam in range(3*len(atoms)):
        G1_ao = np.einsum('uvls,ls->uv', dERI[lam], P)
        K1_ao = np.einsum('ulvs,ls->uv', dERI[lam], P)
        G1_ao = 2*G1_ao - K1_ao
        
        F1_ao = dH[lam] + G1_ao
        F1_mo = C.T @ F1_ao @ C
        S1_mo = C.T @ dS[lam] @ C
        
        RHS = np.zeros(dim)
        for idx in range(dim):
            a, i = map_ai[idx]
            a_mo = n_occ + a
            RHS[idx] = F1_mo[a_mo, i] - eps[i] * S1_mo[a_mo, i]
            
        try:
            U_vec = solve(A_mat, RHS)
        except LinAlgError:
            print("Warning: Singular CPHF matrix, using regularization.")
            U_vec = solve(A_mat + 1e-8*np.eye(dim), RHS)
            
        for idx in range(dim):
            a, i = map_ai[idx]
            U[lam, a, i] = U_vec[idx]

    # 6. Force Constant Assembly (Bond Stretch)
    idx_target = 3 + 2 # Atom 1, Z
    
    d2Enuc_target = 0.0
    if len(atoms)==2: 
        d2Enuc_target = 2 * atoms[0].Z * atoms[1].Z / (r_bond**3)

    G2_ao = np.einsum('uvls,ls->uv', d2ERI[idx_target, idx_target], P)
    K2_ao = np.einsum('ulvs,ls->uv', d2ERI[idx_target, idx_target], P)
    G2_ao = 2*G2_ao - K2_ao
    
    W = np.zeros((n_basis, n_basis))
    for i in range(n_occ):
        W += 2 * eps[i] * np.outer(C[:,i], C[:,i])
        
    term_static = d2Enuc_target + np.sum(P * d2H[idx_target, idx_target]) + \
                  0.5 * np.sum(P * G2_ao) - np.sum(W * d2S[idx_target, idx_target])
    
    G1_ao = np.einsum('uvls,ls->uv', dERI[idx_target], P)
    K1_ao = np.einsum('ulvs,ls->uv', dERI[idx_target], P)
    G1_ao = 2*G1_ao - K1_ao
    F1_mo = C.T @ (dH[idx_target] + G1_ao) @ C
    S1_mo = C.T @ dS[idx_target] @ C
    
    relax_sum = 0.0
    for a in range(n_vir):
        for i in range(n_occ):
            a_mo = n_occ + a
            rhs_val = F1_mo[a_mo, i] - eps[i] * S1_mo[a_mo, i]
            relax_sum += U[idx_target, a, i] * rhs_val
            
    term_relax = 4 * relax_sum 
    
    k_tot = term_static + term_relax
    print(f"  Analytic Force Constant (k_zz): {k_tot:.6f} Eh/bohr^2")
    print(f"    Static: {term_static:.6f}")
    print(f"    Relax : {term_relax:.6f}")

    if verify:
        _verify_hessian_values(atoms, d2S, d2ERI)

def _verify_hessian_values(atoms, d2S, d2ERI):
    print("\n--- Numerical Verification (Integral Values) ---")
    h = 0.001
    
    def get_ints(shift_list):
        origins = [at.coords.copy() for at in atoms]
        basis_orgs = [[b.origin.copy() for b in at.basis] for at in atoms]
        
        for (idx, ax, val) in shift_list:
            atoms[idx].coords[ax] += val
            for b in atoms[idx].basis: b.origin[ax] += val
            
        eng = IntegralEngine(atoms, [b for a in atoms for b in a.basis])
        S, _, _, ERI = eng.compute_integrals()
        
        for i, at in enumerate(atoms):
            at.coords = origins[i]
            for j, b in enumerate(at.basis):
                b.origin = basis_orgs[i][j]
        return S, ERI

    # Mixed deriv check: H1_z, H2_z
    idx1, ax1 = 0, 2
    idx2, ax2 = 1, 2
    
    S_pp, ERI_pp = get_ints([(idx1, ax1, h), (idx2, ax2, h)])
    S_pm, ERI_pm = get_ints([(idx1, ax1, h), (idx2, ax2, -h)])
    S_mp, ERI_mp = get_ints([(idx1, ax1, -h), (idx2, ax2, h)])
    S_mm, ERI_mm = get_ints([(idx1, ax1, -h), (idx2, ax2, -h)])
    
    num_S = (S_pp - S_pm - S_mp + S_mm) / (4*h**2)
    num_ERI = (ERI_pp - ERI_pm - ERI_mp + ERI_mm) / (4*h**2)
    
    anl_S = d2S[2, 5]
    anl_ERI = d2ERI[2, 5]
    
    print(f"Overlap <s1|s2>: Anl={anl_S[0,1]:.8f}, Num={num_S[0,1]:.8f}")
    print(f"ERI (11|22)    : Anl={anl_ERI[0,0,1,1]:.8f}, Num={num_ERI[0,0,1,1]:.8f}")
    
    # Automated assertions
    assert abs(anl_S[0,1] - num_S[0,1]) < 1e-4, "Overlap verification failed!"
    assert abs(anl_ERI[0,0,1,1] - num_ERI[0,0,1,1]) < 1e-4, "ERI verification failed!"
    print("SUCCESS: Analytic 2nd derivatives match numerical results.")

if __name__ == "__main__":
    run_calculation('H2', 1.4, verify=True)

cphf_extended.py

import numpy as np
from scipy.linalg import solve, LinAlgError
import cphf

# =============================================================================
# 1. Extended Engines (Dipole & Gradients)
# =============================================================================

class PropertyIntegralEngine(cphf.IntegralEngine):
    def compute_dipole_integrals(self):
        """Computes <mu|r|nu>. Returns (3, N, N) array."""
        n = self.n_basis
        D_ao = np.zeros((3, n, n)) 
        for i in range(n):
            for j in range(i, n):
                val = self._dipole_element(self.basis[i], self.basis[j])
                for dim in range(3):
                    D_ao[dim, i, j] = val[dim]
                    if i != j: D_ao[dim, j, i] = val[dim]
        return D_ao

    def _dipole_element(self, bi, bj):
        vals = np.zeros(3)
        vec_AB = bi.origin - bj.origin
        dist_sq = np.dot(vec_AB, vec_AB)
        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b
                rp = (a * bi.origin + b * bj.origin) / p
                ov = (np.pi / p)**1.5 * np.exp(-(a * b / p) * dist_sq)
                vals += ca * cb * ov * rp
        return vals

class PropertyGradient(cphf.AnalyticGradient):
    def compute_dipole_gradients(self):
        """Computes d/dR <mu|r|nu>."""
        n_coords = 3 * self.n_atoms
        dD = np.zeros((n_coords, 3, self.n_basis, self.n_basis)) 
        for i in range(self.n_basis):
            for j in range(i, self.n_basis):
                bi, bj = self.basis[i], self.basis[j]
                idx_A, idx_B = self.basis_atom_map[i], self.basis_atom_map[j]
                grads = self._grad_dipole(bi, bj, idx_A, idx_B)
                for at_idx, g_vecs in grads.items():
                    for dim_d in range(3):
                        for dim_r in range(3):
                            row = at_idx*3 + dim_r
                            val = g_vecs[dim_d][dim_r]
                            dD[row, dim_d, i, j] += val
                            if i != j: dD[row, dim_d, j, i] += val
        return dD

    def _grad_dipole(self, bi, bj, idx_A, idx_B):
        res = {idx_A: [np.zeros(3) for _ in range(3)], idx_B: [np.zeros(3) for _ in range(3)]}
        A, B = bi.origin, bj.origin
        vec_AB = A - B
        R2 = np.dot(vec_AB, vec_AB)
        for ia, ca in enumerate(bi.norm_coefs):
            for jb, cb in enumerate(bj.norm_coefs):
                a, b = bi.exps[ia], bj.exps[jb]
                p = a + b; xi = a*b/p; rp = (a*A + b*B)/p
                S_prim = (np.pi/p)**1.5 * np.exp(-xi * R2)
                dS_dA = -2 * xi * vec_AB * S_prim
                term = ca * cb
                for k in range(3):
                    grad_A = dS_dA * rp[k]; grad_A[k] += S_prim * (a/p)
                    grad_B = -dS_dA * rp[k]; grad_B[k] += S_prim * (b/p)
                    res[idx_A][k] += term * grad_A; res[idx_B][k] += term * grad_B
        return res

# =============================================================================
# 2. Main Analysis Logic
# =============================================================================

def get_scf_and_properties(atoms, basis_flat, n_occ):
    """Runs SCF and computes basic properties needed for CPHF."""
    eng = PropertyIntegralEngine(atoms, basis_flat)
    S, T, V, ERI = eng.compute_integrals()
    D_ao = eng.compute_dipole_integrals()
    
    Hcore = T + V
    evals, evecs = np.linalg.eigh(S)
    X = evecs @ np.diag(1.0/np.sqrt(np.clip(evals, 1e-8, None))) @ evecs.T
    
    P = np.zeros((len(basis_flat), len(basis_flat)))
    C = np.zeros_like(P)
    eps = np.zeros(len(basis_flat))
    
    for _ in range(50):
        G = np.einsum('uvls,ls->uv', ERI, P)
        K = np.einsum('ulvs,ls->uv', ERI, P)
        F = Hcore + 2*G - K
        eps, C_p = np.linalg.eigh(X.T @ F @ X)
        C = X @ C_p
        P = C[:, :n_occ] @ C[:, :n_occ].T
        
    return C, eps, P, S, ERI, D_ao

def build_hessian_A(C, eps, ERI, n_occ):
    """Builds the Electronic Hessian (A matrix) for CPHF."""
    n_basis = C.shape[0]
    n_vir = n_basis - n_occ
    
    # Transform ERI to MO
    tmp = np.einsum('uvls, up->pvls', ERI, C)
    tmp = np.einsum('pvls, vq->pqls', tmp, C)
    tmp = np.einsum('pqls, lr->pqrs', tmp, C)
    I_mo = np.einsum('pqrs, st->pqrt', tmp, C)
    
    dim = n_occ * n_vir
    A_mat = np.zeros((dim, dim))
    map_ai = []
    for a in range(n_vir):
        for i in range(n_occ):
            map_ai.append((a, i))
            
    for idx1 in range(dim):
        a, i = map_ai[idx1]
        a_mo, i_mo = n_occ + a, i
        for idx2 in range(dim):
            b, j = map_ai[idx2]
            b_mo, j_mo = n_occ + b, j
            val = (eps[a_mo] - eps[i_mo]) if idx1 == idx2 else 0.0
            val += 4 * I_mo[a_mo, i_mo, b_mo, j_mo]
            val -= I_mo[a_mo, b_mo, j_mo, i_mo]
            val -= I_mo[a_mo, j_mo, b_mo, i_mo]
            A_mat[idx1, idx2] = val
            
    return A_mat, map_ai, I_mo

def solve_dipole_derivative(atoms, basis_flat, n_occ, target_idx, target_axis):
    """Calculates Analytic Dipole Derivative (IR)."""
    C, eps, P, S, ERI, D_ao = get_scf_and_properties(atoms, basis_flat, n_occ)
    A_mat, map_ai, _ = build_hessian_A(C, eps, ERI, n_occ)
    
    # Geometric Derivatives
    ag = PropertyGradient(atoms, basis_flat)
    dS, dH, dERI = ag.run() # Now using FIXED cphf.py
    dD_ao = ag.compute_dipole_gradients()
    
    pert_idx = target_idx * 3 + target_axis
    
    # RHS for Geometric CPHF
    G1 = 2*np.einsum('uvls,ls->uv', dERI[pert_idx], P) - np.einsum('ulvs,ls->uv', dERI[pert_idx], P)
    F1_mo = C.T @ (dH[pert_idx] + G1) @ C
    S1_mo = C.T @ dS[pert_idx] @ C
    
    dim = len(map_ai)
    RHS = np.zeros(dim)
    for idx in range(dim):
        a, i = map_ai[idx]
        RHS[idx] = -(F1_mo[n_occ+a, i] - eps[i] * S1_mo[n_occ+a, i])
        
    try:
        U_geo = solve(A_mat, RHS)
    except LinAlgError:
        U_geo = solve(A_mat + 1e-8*np.eye(dim), RHS)
        
    # Dipole Derivative Assembly
    dmu_dR = np.zeros(3)
    dmu_dR[target_axis] += atoms[target_idx].Z # Nuclear
    
    D_mo = [C.T @ D_ao[k] @ C for k in range(3)]
    
    for k in range(3):
        # 1. Hellmann-Feynman
        term = np.sum(P * dD_ao[pert_idx, k])
        
        # 2. Relaxation (Mixing)
        for idx in range(dim):
            a, i = map_ai[idx]
            term += 2.0 * U_geo[idx] * D_mo[k][n_occ+a, i]
            
        # 3. Renormalization (Occ-Occ) 
        for i in range(n_occ):
            for j in range(n_occ):
                term += (-S1_mo[i, j]) * D_mo[k][j, i]
                
        dmu_dR[k] -= 2.0 * term # Elec coeff
        
    return dmu_dR

def solve_polarizability(atoms, basis_flat, n_occ):
    """Calculates Analytic Static Polarizability Tensor (alpha)."""
    C, eps, P, S, ERI, D_ao = get_scf_and_properties(atoms, basis_flat, n_occ)
    A_mat, map_ai, _ = build_hessian_A(C, eps, ERI, n_occ)
    
    D_mo = [C.T @ D_ao[k] @ C for k in range(3)]
    n_vir = len(basis_flat) - n_occ
    dim = len(map_ai)
    
    # Solve CPHF for Electric Field (x, y, z)
    # H^(1) = r (Dipole operator)
    # RHS_ai = - <a|r|i>
    
    alpha = np.zeros((3, 3))
    
    for comp_field in range(3):
        RHS_elec = np.zeros(dim)
        for idx in range(dim):
            a, i = map_ai[idx]
            RHS_elec[idx] = - D_mo[comp_field][n_occ+a, i]
            
        try:
            U_elec = solve(A_mat, RHS_elec)
        except LinAlgError:
            U_elec = solve(A_mat + 1e-8*np.eye(dim), RHS_elec)
            
        # Alpha_uv = - Tr(mu_u * P^(v))
        #          = - 4 * sum_ai U^v_ai * <a|mu_u|i>
        for comp_resp in range(3):
            val = 0.0
            for idx in range(dim):
                a, i = map_ai[idx]
                val += U_elec[idx] * D_mo[comp_resp][n_occ+a, i]
            alpha[comp_resp, comp_field] = -4.0 * val
            
    return alpha

# =============================================================================
# 3. Full Spectra Driver
# =============================================================================

def run_spectroscopy_analysis(mol_name, R0):
    print(f"\n{'='*60}")
    print(f" SPECTROSCOPIC ANALYSIS: {mol_name} (R={R0:.4f})")
    print(f"{'='*60}")
    
    # 1. Define Molecule System
    def get_system(r):
        if mol_name == 'H2':
            a1 = cphf.Atom('H', [0,0,0], 1)
            a2 = cphf.Atom('H', [0,0,r], 1)
            exps = [3.42525, 0.623913, 0.168855]
            coefs = [0.154329, 0.535328, 0.444635]
            a1.basis.append(cphf.BasisFunction(a1.coords, (exps, coefs)))
            a2.basis.append(cphf.BasisFunction(a2.coords, (exps, coefs)))
            return [a1, a2], 1, 2 # Moving Atom Index, Axis (Z)
        elif mol_name == 'LiH':
            a1 = cphf.Atom('Li', [0,0,0], 3)
            a2 = cphf.Atom('H', [0,0,r], 1)
            
            # Li 1s (STO-3G)
            exps_li_1s = [16.119575, 2.9362007, 0.7946505]
            coefs_li_1s = [0.15432897, 0.53532814, 0.44463454]
            a1.basis.append(cphf.BasisFunction(a1.coords, (exps_li_1s, coefs_li_1s)))
            
            # Li 2s (STO-3G)
            exps_li_2s = [0.6362897, 0.1478601, 0.0480887]
            coefs_li_2s = [-0.09996723, 0.39951283, 0.70011547]
            a1.basis.append(cphf.BasisFunction(a1.coords, (exps_li_2s, coefs_li_2s)))
            
            # H 1s (STO-3G)
            exps_h = [3.42525, 0.623913, 0.168855]
            coefs_h = [0.154329, 0.535328, 0.444635]
            a2.basis.append(cphf.BasisFunction(a2.coords, (exps_h, coefs_h)))
            
            # System: Li(1s,2s) + H(1s) = 3 basis functions
            # Electrons: 4 => N_occ=2, N_vir=1
            return [a1, a2], 1, 2

    atoms, idx, axis = get_system(R0)
    basis = [b for a in atoms for b in a.basis]
    n_occ = sum(a.Z for a in atoms)//2
    
    # 2. IR: Dipole Derivative (Analytic)
    dmu_dR = solve_dipole_derivative(atoms, basis, n_occ, idx, axis)
    ir_intensity = np.linalg.norm(dmu_dR)**2 # Simple squared magnitude
    
    print(f"\n[IR Spectrum]")
    print(f"  Dipole Deriv (vector): {dmu_dR}")
    print(f"  IR Intensity (|dmu/dR|^2): {ir_intensity:.6e}")
    if ir_intensity < 1e-6:
        print("  >> IR INACTIVE (Forbidden)")
    else:
        print("  >> IR ACTIVE")

    # 3. Raman: Polarizability Derivative (Analytic Alpha + Finite Diff)
    # Step A: Analytic Alpha at R0
    alpha_0 = solve_polarizability(atoms, basis, n_occ)
    
    # Step B: Analytic Alpha at R0 + h
    h = 0.001
    atoms_p, _, _ = get_system(R0 + h)
    basis_p = [b for a in atoms_p for b in a.basis]
    alpha_p = solve_polarizability(atoms_p, basis_p, n_occ)
    
    # Derivative
    d_alpha_dR = (alpha_p - alpha_0) / h
    
    # Raman Activity (Simplified for diatomic stretch)
    # Usually involves mean polarizability deriv (a') and anisotropy deriv (gamma')
    # Here we just show the tensor derivative
    print(f"\n[Raman Spectrum]")
    print(f"  Static Polarizability (R={R0}):\n{alpha_0}")
    print(f"  d(Alpha)/dR (Finite Diff of Analytic Alpha):\n{d_alpha_dR}")
    
    # For linear molecule along Z:
    # Parallel component (zz) and Perpendicular (xx, yy)
    da_par = d_alpha_dR[2, 2]
    da_perp = d_alpha_dR[0, 0]
    
    raman_activity = da_par**2 + 2 * da_perp**2 # Approx isotropic measure
    print(f"  Raman Activity (~ da_par^2 + 2*da_perp^2): {raman_activity:.6e}")
    
    if raman_activity < 1e-6:
        print("  >> RAMAN INACTIVE")
    else:
        print("  >> RAMAN ACTIVE")

if __name__ == "__main__":
    run_spectroscopy_analysis('H2', 1.4)
    run_spectroscopy_analysis('LiH', 3.015)