Coupled-Perturbed Hartree-Fock (CPHF) 法の数理と実装：分極率・双極子モーメント微分の厳密な導出

最終更新：2025-12-30

注意: この記事はAIによって自動生成されたものです。正確な情報については、Szabo & Ostlundの「Modern Quantum Chemistry」や原著論文をご確認ください。

序論：動的な電子状態への応答#

Hartree-Fock（HF）法は、全エネルギーを分子軌道（MO）係数の関数として最小化する静的な手法である。しかし、我々が実験的に観測する多くの物理量は、分子の「変化」に対する応答である。

振動数（Hessian）: 核の微小変位に対するエネルギーの応答。
分極率（Polarizability）: 外部電場の微小変化に対する誘起双極子モーメントの応答。
赤外（IR）強度: 核の微小変位に対する双極子モーメントの応答。

これらの量を正確に評価するためには、単に「未摂動の電子密度」を用いるだけでは不十分であり、摂動（核変位や電場）によってMO係数自体がどのように弛緩（リラックス）するかを決定しなければならない。この応答を線形方程式として定式化したものが、Coupled-Perturbed Hartree-Fock (CPHF) 法である。

本稿では、CPHFの数理的導出を極めて詳細に展開し、実際のプログラム実装に直結するアルゴリズムを解説する。

1. 数理的背景：摂動下での変分条件#

1.1 MO係数の微分と $U$ 行列の定義#

パラメータ $\lambda$ （電場や核座標）を摂動とすると、摂動下のMO係数 $\mathbf{C}(\lambda)$ は、未摂動のMO係数 $\mathbf{C}^{(0)}$ を用いて次のように展開できる。

$\mathbf{C}^\lambda = \mathbf{C} \mathbf{U}^\lambda \tag{1}$

ここで $\mathbf{U}^\lambda$ は軌道間の混合を記述する行列である。 MOの正規直交化条件 $\mathbf{C}^\dagger \mathbf{S} \mathbf{C} = \mathbf{I}$ を微分すると、次式の制約が得られる。

$\mathbf{U}^{\lambda \dagger} + \mathbf{U}^\lambda + \mathbf{S}^{\lambda, \text{MO}} = 0 \tag{2}$

ここで $S_{pq}^{\lambda, \text{MO}}$ は重なり行列のMO基底での微分である。対角成分 $U_{ii}^\lambda$ は $- \frac{1}{2} S_{ii}^{\lambda, \text{MO}}$ と即座に求まるが、非対角成分、特に占有軌道 $i$ と仮想軌道 $a$ の間の $U_{ai}^\lambda$ はFock行列の応答を介して複雑に結合する。

1.2 Fock方程式の一次微分#

収束したSCF状態では $F_{ai}^{\text{MO}} = 0$ である。摂動下でもこの条件が維持されるためには、その全微分がゼロでなければならない。

$\frac{d}{d\lambda} F_{ai}^{\text{MO}} = F_{ai}^{\lambda, \text{MO}} + (\epsilon_i - \epsilon_a) U_{ai}^\lambda = 0 \tag{3}$

ここで、Fock行列の微分 $F_{ai}^{\lambda, \text{MO}}$ には、密度行列の微分 $\mathbf{P}^\lambda$ を通じて $U_{ai}^\lambda$ 自身が含まれる。 $\mathbf{P}^\lambda = 2 \sum_i^{occ} (\mathbf{C} \mathbf{U}_i^\lambda \mathbf{C}_i^\dagger + \mathbf{C}_i \mathbf{U}_i^{\lambda \dagger} \mathbf{C}^\dagger)$

これを代入し、整理すると、占有・仮想ブロックに関するCPHF方程式が得られる。

$(\epsilon_a - \epsilon_i) U_{ai}^\lambda + \sum_{bj} [4(ai|bj) - (ab|ij) - (aj|ib)] U_{bj}^\lambda = B_{ai}^\lambda \tag{4}$

ここで $B_{ai}^\lambda$ は $U$ を含まない項（積分微分など）である。

2. 数理的導出：CPHF方程式#

本節では、閉殻系の制限Hartree-Fock (RHF) 法を対象として、パラメータ $\lambda$ に対する応答を導出する。

2.1 未摂動系の条件#

HF方程式を行列形式（Roothaan方程式）で示す： $\mathbf{F}\mathbf{C} = \mathbf{S}\mathbf{C}\boldsymbol{\epsilon}$ ここで $\mathbf{F}$ はFock行列、 $\mathbf{C}$ はMO係数、 $\mathbf{S}$ は重なり行列、 $\boldsymbol{\epsilon}$ は軌道エネルギー対角行列である。また、MOの正規直交化条件は以下である： $\mathbf{C}^\dagger \mathbf{S} \mathbf{C} = \mathbf{1}$

2.2 摂動下での変分条件#

パラメータ $\lambda$ が微小変化したとき、各行列も $\lambda$ に依存して変化する。 $\frac{d}{d\lambda} (\mathbf{C}^\dagger \mathbf{S} \mathbf{C}) = \mathbf{C}^{\lambda \dagger} \mathbf{S} \mathbf{C} + \mathbf{C}^\dagger \mathbf{S}^\lambda \mathbf{C} + \mathbf{C}^\dagger \mathbf{S} \mathbf{C}^\lambda = 0$ ここで肩付きの $\lambda$ は $\lambda$ に関する微分を表す。

MO係数の変化 $\mathbf{C}^\lambda$ を、元の係数 $\mathbf{C}$ の線形結合で表す（ユニタリ変換の拡張）： $\mathbf{C}^\lambda = \mathbf{C} \mathbf{U}^\lambda$ この $\mathbf{U}^\lambda$ 行列を決定することがCPHFの目的である。正規直交化条件の微分にこれを代入すると、 $\mathbf{U}^{\lambda \dagger} + \mathbf{U}^\lambda + \mathbf{S}^{\lambda, \text{MO}} = 0$ ここで $\mathbf{S}^{\lambda, \text{MO}} = \mathbf{C}^\dagger \mathbf{S}^\lambda \mathbf{C}$ である。対角成分については $2 U_{ii}^\lambda + S_{ii}^{\lambda, \text{MO}} = 0$ より直ちに求まる。

2.3 Fock行列の応答#

Fock行列の微分 $\mathbf{F}^\lambda$ は、以下の二つの寄与からなる：

Core寄与: 1電子積分や基底関数微分の直接的変化。
密度変化寄与: $\mathbf{C}^\lambda$ を通じた2電子積分項の変化。

MO基底でのFock行列の微分を考えると、正準MO条件下（ $\mathbf{F}_{MO}$ が対角行列）において、 $\mathbf{F}\mathbf{C} = \mathbf{S}\mathbf{C}\boldsymbol{\epsilon}$ の微分から次の方程式が得られる： $F_{pq}^{\lambda, \text{MO}} - S_{pq}^{\lambda, \text{MO}} \epsilon_q - \delta_{pq} \epsilon_p^\lambda = (\epsilon_q - \epsilon_p) U_{pq}^\lambda$

ここから、占有軌道 $i$ と仮想軌道 $a$ の間の成分 $U_{ai}^\lambda$ についての結合連立方程式が導かれる。

2.4 CPHF連立一次方程式#

最終的に、 $U_{ai}^\lambda$ は以下の形式の線形方程式の解として得られる： $(\epsilon_a - \epsilon_i) U_{ai}^\lambda + \sum_{bj} A_{ai, bj} U_{bj}^\lambda = B_{ai}^\lambda$ ここで、行列 $\mathbf{A}$ は軌道間相互作用（Hessianの一部）を、 $\mathbf{B}$ は摂動の直接項を表す。 $A_{ai, bj} = 4(ai|bj) - (ab|ij) - (aj|ib)$ $B_{ai}^\lambda$ は、MO係数の微分を含まない項（Fock行列の直接微分など）を集約したものである。

3. 物理量への応用：分極率と双極子モーメント微分#

3.1 外部電場に対する応答（分極率）#

外部電場 $\mathcal{E}_z$ に対する応答を考える場合、摂動項は双極子行列 $\mathbf{D}_z$ となる（ $H^\lambda = -D_z$ ）。このとき、重なり行列の微分 $S^\lambda$ はゼロとなるため、式(2)より $U_{ii} = 0$ 、かつ $U_{ij} = -U_{ji}$ となる。

分極率 $\alpha_{zz}$ は次のように定義される。 $\alpha_{zz} = \frac{d\mu_z}{d\mathcal{E}_z} = \sum_{\mu\nu} P_{\mu\nu} \frac{dD_{z, \mu\nu}}{d\mathcal{E}_z} + \sum_{\mu\nu} \frac{dP_{\mu\nu}}{d\mathcal{E}_z} D_{z, \mu\nu}$ ここで第2項が電子状態の弛緩（CPHFの解 $U^\lambda$ ）による寄与であり、これを含めることで初めて精度の高い分極率が得られる。

3.2 核変位に対する応答（双極子モーメント微分）#

赤外強度の計算に必要な $d\boldsymbol{\mu}/dR_A$ においても同様である。核が動くと基底関数が動くため、積分の微分 $\langle \chi_\mu^\lambda | r | \chi_\nu \rangle$ に加え、電子密度の変化（CPHF）が寄与する。これは、外部電場をかけずとも、核の動きによって「電子が引きずられる」効果を厳密に計算していることに相当する。

4. 実装詳細：反復法アルゴリズム#

CPHF方程式 $( \mathbf{A} - \Delta\boldsymbol{\epsilon} ) \mathbf{U} = \mathbf{B}$ を解くためには、行列 $\mathbf{A}$ を直接構成するのではなく、応答 Fock 行列を AO 基底で逐次生成する反復法が効率的である。

以下に、CPHFによる分極率計算の完全な実装ロジックを示す。

import numpy as np
import copy

def run_cprhf(molecular_orbitals, orbital_energy, dipole_z, eri_ao, nelec, n_basisset, max_iter=50, tol=1e-10):
    """
    Coupled-Perturbed Restricted Hartree-Fock (CPRHF) Implementation
    for calculating polarizability alpha_zz.
    """
    n_occ = int(nelec / 2)
    n_virt = n_basisset - n_occ
    
    # MO indices
    occ = slice(0, n_occ)
    virt = slice(n_occ, n_basisset)

    # 1. Transform Dipole Matrix to MO basis (B_ai term)
    # Dipole matrix in MO: <phi_p | z | phi_q>
    dip_mo = np.dot(molecular_orbitals.T, np.dot(dipole_z, molecular_orbitals))
    
    # RHS: B_ai = - <phi_a | z | phi_i>
    # (assuming S_lambda = 0 for electric field perturbation)
    B = -1.0 * dip_mo[virt, occ]

    # 2. Initial Guess: Uncoupled Approximation (U = B / delta_epsilon)
    epsilon_occ = orbital_energy[occ]
    epsilon_virt = orbital_energy[virt]
    denom = epsilon_occ[None, :] - epsilon_virt[:, None] # eps_i - eps_a
    
    U_mo = B / denom
    U_mo_old = copy.deepcopy(U_mo)

    print("Starting CPRHF iterations...")
    
    for itr in range(max_iter):
        # 3. Transform U_mo back to AO basis to construct response density
        # P_resp = sum_{ai} (C_a U_ai C_i^T + C_i U_ai^T C_a^T) * 2.0
        # This is the "Change in Density Matrix" due to the field
        C_occ = molecular_orbitals[:, occ]
        C_virt = molecular_orbitals[:, virt]
        
        # AO response density matrix
        P_resp = np.dot(np.dot(C_virt, U_mo), C_occ.T)
        P_resp = (P_resp + P_resp.T) * 2.0 # Symmetry and occupancy factor

        # 4. Construct Response J and K matrices in AO basis
        # This is the "Coupled" part: how density change affects the potential
        J_resp, K_resp = construct_jk_from_eri(eri_ao, P_resp)
        
        # 5. Transform Response Fock to MO basis
        F_resp_ao = J_resp - 0.5 * K_resp
        F_resp_mo = np.dot(C_virt.T, np.dot(F_resp_ao, C_occ))

        # 6. Update U_mo using the CPHF Equation
        # U_ai = (B_ai - F_resp_mo_ai) / (eps_i - eps_a)
        U_mo = (B + F_resp_mo) / denom

        # 7. Convergence Check
        residual = np.linalg.norm(U_mo - U_mo_old)
        print(f"Iteration {itr}: Residual = {residual:.12e}")
        
        if residual < tol:
            print("CPRHF Converged.")
            break
        U_mo_old = copy.deepcopy(U_mo)
    else:
        print("CPRHF failed to converge.")

    # 8. Calculate Polarizability alpha_zz
    # alpha = 2 * sum_{ai} U_ai * <phi_a | z | phi_i> * 2.0 (spin)
    # Or in AO basis using the response density:
    alpha_zz = -2.0 * np.sum(P_resp * dipole_z)
    
    return alpha_zz, P_resp

def construct_jk_from_eri(eri, P):
    # Standard JK construction from 4-center ERI and density P
    n = P.shape[0]
    J = np.einsum('pqrs,rs->pq', eri, P)
    K = np.einsum('prqs,rs->pq', eri, P)
    return J, K

数理的な行間の補足#

分母 $D_{ai}$ : $U_{ai}$ は軌道エネルギー差 $(\epsilon_i - \epsilon_a)$ に反比例する。これは、HOMO-LUMOギャップが小さい系ほど電場応答（分極率）が大きくなるという物理的実態を数理的に示している。

結合（Coupled）項: J_resp - 0.5 * K_resp の項が存在しない場合、これはUncoupled Hartree-Fock (UCHF) と呼ばれる。

CPHFは、電子密度の偏りが自分自身の作るポテンシャルを変化させ、それがさらに電子密度の偏りを誘強するという自己矛盾のない応答を記述している。

補足#

「風」による最初の刺激（右辺の作成）

物理的なイメージ：「撮影中に横から強い風（電場）が吹いてきた」

# 1. Transform Dipole Matrix to MO basis (B_ai term)
dip_mo = np.dot(molecular_orbitals.T, np.dot(dipole_z, molecular_orbitals))
B = -1.0 * dip_mo[virt, occ]

コードの意味: dipole_z（電場という刺激）を、現在の分子軌道（MO）の言葉に翻訳しています。

数理の行間: B は、電子を「占有軌道から仮想軌道へ押し出す力」の強さを表しています。電場という外圧が、どれくらい電子の座席をかき乱そうとしているか、その第一歩を定義しています。

「とりあえず傾いてみる」（初期推測値）物理的なイメージ：「隣の人のことは気にせず、風の強さ分だけ傾いてみる（Uncoupled）」

# 2. Initial Guess: Uncoupled Approximation
denom = epsilon_occ[None, :] - epsilon_virt[:, None] 
U_mo = B / denom

コードの意味: 刺激 B を、座席の移動のしにくさ denom（エネルギー差）で割り算しています。

数理の行間: ここで一旦 U_mo（電子の傾き具合）が出ますが、これはまだ「隣の人との反発の変化」を無視した、大雑把な予測値です。

「隣の人の変化」を計算する（応答密度）物理的なイメージ：「自分が傾いたことで、電子の密度の濃い場所が変わる」

# 3. Transform U_mo back to AO basis to construct response density
P_resp = np.dot(np.dot(C_virt, U_mo), C_occ.T)
P_resp = (P_resp + P_resp.T) * 2.0

コードの意味: 軌道の傾き U_mo を、「電子密度の変化（P_resp）」という具体的な形に変換しています。

数理の行間: 電場そのものによる変化ではなく、「電場で電子が動いた結果、新しく生じた密度の偏り」を算出しています。

「押し合い」の発生（J/K行列の構成）物理的なイメージ：「隣の人が傾いてきたせいで、反発が強くなって邪魔になる」

# 4. Construct Response J and K matrices in AO basis
J_resp, K_resp = construct_jk_from_eri(eri_ao, P_resp)

# 5. Transform Response Fock to MO basis
F_resp_ao = J_resp - 0.5 * K_resp
F_resp_mo = np.dot(C_virt.T, np.dot(F_resp_ao, C_occ))

コードの意味: 密度の変化 P_resp が原因で発生する、新しい電子間反発（J）と交換相互作用（K）を計算しています。

数理の行間: ここが Coupled（結合）の核心です。電子が動いたことで「景色（ポテンシャル）」が変わり、それが自分たちに跳ね返ってくるプロセスです。F_resp_mo は「自分たちの動きによって生じた、自分たちへの新しい壁」です。

「ちょうどいい姿勢」への修正（反復計算）物理的なイメージ：「風の力」と「隣の人からの反発」のバランスが取れるまで微調整を繰り返す

# 6. Update U_mo using the CPHF Equation
U_mo = (B + F_resp_mo) / denom

コードの意味: 「元の風の力 B」に、「新しく生まれた反発 F_resp_mo」を足して、もう一度姿勢 U_mo を決め直します。

数理の行間: B（外圧）と F_resp_mo（内圧）の合計が、ちょうど U_mo を生み出す力と釣り合うまで、このループを回し続けます（for itr in range(max_iter)）。

最終的な「踏ん張り」の強さを算出（分極率）物理的なイメージ：「最終的にどれくらい電子が偏ったか（傾いたか）を測定する」

# 8. Calculate Polarizability alpha_zz
alpha_zz = -2.0 * np.sum(P_resp * dipole_z)

コードの意味: 全ての微調整が終わった後の「密度の偏り P_resp」と「電場 dipole_z」を掛け合わせます。

数理の行間: 最終的に「どれだけのエネルギーが電場によって得した（あるいは損した）か」という応答の強さを数値化しています。これが分極率（電場に対する敏感さ）です。

5. 歴史的・実利的な意義#

解析的Hessianの実現: CPHFにより、核座標微分 $\mathbf{U}^R$ が得られるようになったことで、調和振動数の計算が飛躍的に高速化した。それ以前は数値微分（構造を少しずつずらして計算）に頼っており、莫大な計算時間を要していた。

赤外強度の予言: 外部電場がない状態での双極子モーメント微分を厳密に扱うことで、実験スペクトルの強度を定量的、あるいは定性的に正しく再現できるようになった。これは分子構造同定における決定打となった。

分極率の精密評価: 電子相関を取り込む前のベースラインとして、HFレベルでの正確な分極率は物性化学の基礎である。

結論#

CPHF法は、変分条件という静的な制約を摂動論的に拡張し、電子状態の「弛緩」を線形代数の問題に落とし込んだ。導出の過程で見られた $U$ 行列の直交性や Fock 行列の応答は、現代の TD-DFT（時間依存密度汎関数法）の数理構造の雛形でもあり、応答理論を理解する上での最良の教科書といえる。

参考文献#

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 (1968).
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. 16, 225 (1979).
A. Szabo and N. S. Ostlund, “Modern Quantum Chemistry”, Dover Publications (1996).
Y. Yamaguchi et al., “Analytic Derivative Methods in Ab Initio Molecular Electronic Structure Theory”, Oxford University Press (1994).
https://github.com/ss0832/HartreeFock_py/blob/master/CPHF/rhf_h2/cprhf_h2.py (CPHF法を使って水素分子の双極子モーメント等を求めるために実装してみた例)