last_modified: 2026-01-23

注記:これまでの s-s, s-p, s-d の流れを踏まえ、今回は化学結合において最も基本的かつ重要な p-p 重なり（$\sigma$ 結合および $\pi$ 結合） の導出を行います。ここでも「正規化係数」と「未正規化積分」を厳密に分離するスタイルを維持します。 注記2: 生成AIへの指示で自動生成し、検算をしておりますが、内容が誤っている可能性があります。ご了承ください。Gemini 3.0 Pro及びGPT-5.2の性能を確かめる目的も兼ねています。

1. 数理的定式化：p-p 重なりの一元化公式#

中心 $A$ の短縮p軌道 $\Psi_{A, \gamma}$ （成分 $\gamma \in \{x,y,z\}$）と、中心 $B$ の短縮p軌道 $\Psi_{B, \delta}$ （成分 $\delta \in \{x,y,z\}$）の重なり積分 $S_{\gamma\delta}$ を導出します。

公式

各項の定義正規化済み係数 $D$:$D_{A,\mu}(\gamma) = d_{A,\mu} N_{A,\mu}^{cont} N_{prim}(\alpha_\mu, \mathbf{l}_\gamma)$

ここでは p軌道の正規化定数 $N_{prim}$（$l=1$）を使用します。

基本ガウス積分 $S_{ss}^{(\mu\nu, 0)}$:未正規化の s-s 重なり積分（前項までと同様）。

$S_{ss}^{(\mu\nu, 0)} = (\frac{\pi}{p})^{3/2} e^{-\frac{\alpha\beta}{p} |\mathbf{R}_{AB}|^2}, \quad p = \alpha_\mu + \beta_\nu$

幾何学的因子（シフト項）:ガウス積の中心座標 $\mathbf{P} = \frac{\alpha \mathbf{A} + \beta \mathbf{B}}{\alpha + \beta}$ を用いた相対座標項。

$P_\gamma - A_\gamma = \frac{\beta (B_\gamma - A_\gamma)}{p}$ $P_\delta - B_\delta = \frac{\alpha (A_\delta - B_\delta)}{p} = -\frac{\alpha (B_\delta - A_\delta)}{p}$

物理的意味: この2つの項の積が、軌道の「向き」による符号（結合性 vs 反結合性）を決定します。

等方的な広がり項:$\frac{\delta_{\gamma\delta}}{2p}$: 軌道の広がり（分散）に由来する重なり。$\gamma = \delta$ （例：$p_x - p_x$）の時のみ寄与します。これが $\pi$ 結合の主要因となります。

2. Pythonによる実装#

s-d 重なりで整備したユーティリティ関数を再利用し、p-p 重なり（3成分 $\times$ 3成分 = 9要素）を計算するクラスを作成します。

import numpy as np

# --- 1. 共通数学ユーティリティ (前回のコードから継承・再掲) ---

def double_factorial(n):
    if n <= 1: return 1.0
    res = 1.0
    for i in range(n, 0, -2): res *= i
    return res

def get_primitive_norm(zeta: float, l: int, m: int, n: int) -> float:
    """プリミティブ正規化定数 (厳密解)"""
    L_sum = l + m + n
    prefactor = (2.0 * zeta / np.pi)**0.75
    numerator = (8.0 * zeta)**L_sum
    denominator = double_factorial(2*l - 1) * double_factorial(2*m - 1) * double_factorial(2*n - 1)
    return prefactor * np.sqrt(numerator / denominator)

def calculate_unnormalized_primitive_overlap_generic(alpha, beta, lmn_A, lmn_B, R_vec):
    """
    正規化計算(R=0)のための汎用プリミティブ重なり計算。
    今回は p-p (l=1) の正規化に必要。
    """
    p = alpha + beta
    # R_vec must be 0 for normalization check context usually
    # But full generic implementation requires Obara-Saika recurrence.
    # ここでは「同一中心(R=0)かつ同一成分」の正規化用ケースのみ簡易実装
    if np.linalg.norm(R_vec) < 1e-9 and lmn_A == lmn_B:
        lx, ly, lz = lmn_A
        # int x^2lx exp(-p x^2) -> (2lx-1)!! / (2p)^lx * sqrt(pi/p)
        # Total integral = Ix * Iy * Iz
        def i_1d(k): return double_factorial(2*k - 1) / ((2.0*p)**k) * np.sqrt(np.pi/p)
        return i_1d(lx) * i_1d(ly) * i_1d(lz)
    return 0.0 # Fallback for non-implemented cases in this snippet

def get_contraction_normalizer(exponents, coeffs, lmn):
    """短縮係数の正規化"""
    n = len(exponents)
    s_sum = 0.0
    for i in range(n):
        for j in range(n):
            Ni = get_primitive_norm(exponents[i], *lmn)
            Nj = get_primitive_norm(exponents[j], *lmn)
            S_prim = calculate_unnormalized_primitive_overlap_generic(
                exponents[i], exponents[j], lmn, lmn, np.zeros(3)
            )
            s_sum += coeffs[i] * coeffs[j] * Ni * Nj * S_prim
    return 1.0 / np.sqrt(s_sum)

# --- 2. p-p 重なり積分計算クラス ---

def calculate_overlap_pp_rigorous(
    center_A: np.ndarray, exps_A: list, coeffs_A: list, # p-orbital A
    center_B: np.ndarray, exps_B: list, coeffs_B: list  # p-orbital B
) -> dict:
    """
    p軌道A (px, py, pz) と p軌道B (px, py, pz) の重なり積分行列を計算する。
    
    Returns:
        dict: keys 'xx', 'xy', 'xz', 'yx', ... (Aの成分 + Bの成分)
    """
    R_vec_AB = center_A - center_B # A - B
    # R_vec_BA = center_B - center_A # B - A
    
    # 成分定義
    p_comps = ['x', 'y', 'z']
    lmn_map = {'x':(1,0,0), 'y':(0,1,0), 'z':(0,0,1)}
    
    # 1. 正規化 (A, Bともにp軌道なので、成分ごとの正規化定数は等しいが、念のため計算)
    # p軌道は球対称性により x,y,z で正規化定数は同じ。代表してxで計算。
    N_cont_A = get_contraction_normalizer(exps_A, coeffs_A, (1,0,0))
    N_cont_B = get_contraction_normalizer(exps_B, coeffs_B, (1,0,0))
    
    results = {}
    
    # 2. 二重ループ
    for i, alpha in enumerate(exps_A):
        # Aの係数項 (N_primは成分によらず同じだが、論理的整合性のため関数を通す)
        # N_prim(alpha, 1,0,0) = N_prim(alpha, 0,1,0) ...
        Da_base = coeffs_A[i] * N_cont_A * get_primitive_norm(alpha, 1,0,0)
        
        for j, beta in enumerate(exps_B):
            Db_base = coeffs_B[j] * N_cont_B * get_primitive_norm(beta, 1,0,0)
            
            # 共通項
            p = alpha + beta
            sss_0 = (np.pi / p)**1.5 * np.exp( - (alpha * beta / p) * np.dot(R_vec_AB, R_vec_AB) )
            inv_2p = 1.0 / (2.0 * p)
            
            # ガウス積中心 P の相対座標計算
            # P = (alpha*A + beta*B) / p
            # P - A = beta(B - A) / p = - beta * (A - B) / p
            PA_vec = - (beta / p) * R_vec_AB
            # P - B = alpha(A - B) / p
            PB_vec = (alpha / p) * R_vec_AB
            
            # 9成分の計算
            for comp_a in p_comps: # Aの成分 (gamma)
                for comp_b in p_comps: # Bの成分 (delta)
                    key = comp_a + comp_b
                    
                    # インデックス変換 (x->0, y->1, z->2)
                    idx_a = {'x':0, 'y':1, 'z':2}[comp_a]
                    idx_b = {'x':0, 'y':1, 'z':2}[comp_b]
                    
                    # 数式: (P_gamma - A_gamma)*(P_delta - B_delta) + delta_gd / 2p
                    term_geom = PA_vec[idx_a] * PB_vec[idx_b]
                    
                    if comp_a == comp_b:
                        term_geom += inv_2p
                        
                    # 加算
                    if key not in results: results[key] = 0.0
                    results[key] += Da_base * Db_base * sss_0 * term_geom

    return results

# --- 検証 ---
if __name__ == "__main__":
    # 配置: z軸上の結合 (A=Origin, B=on Z axis)
    pos_A = np.array([0.0, 0.0, 0.0])
    pos_B = np.array([0.0, 0.0, 3.0]) # 距離 3.0 Bohr
    
    # パラメータ (ダミー: 指数1.0)
    exps = [1.0]
    coeffs = [1.0]
    
    overlaps = calculate_overlap_pp_rigorous(pos_A, exps, coeffs, pos_B, exps, coeffs)
    
    print(f"p-p Overlap (Z-axis separation, R=3.0)")
    
    # 行列形式で表示
    mat = np.zeros((3,3))
    keys = ['x', 'y', 'z']
    print("      px_B      py_B      pz_B")
    for i, ca in enumerate(keys):
        row_str = f"p{ca}_A  "
        for j, cb in enumerate(keys):
            val = overlaps[ca+cb]
            mat[i,j] = val
            row_str += f"{val: .5f}  "
        print(row_str)

    print("\n--- Physical Interpretation ---")
    print(f"Sigma Bond (pz-pz): {overlaps['zz']:.5f}")
    print("  -> (PA_z * PB_z) is dominant. Signs: PA_z > 0, PB_z < 0 (pointing inward). Product < 0.")
    print("  -> Plus 1/2p term. Net result depends on distance.")
    
    print(f"Pi Bond (px-px):    {overlaps['xx']:.5f}")
    print("  -> PA_x = 0, PB_x = 0. Geometric term is 0.")
    print("  -> Only 1/2p term contributes. Pure 'sideways' overlap.")
    
    print(f"Orthogonal (px-py): {overlaps['xy']:.5f}")
    print("  -> 0.0 expected by symmetry.")

3. 結果の物理的解釈#

出力結果から、化学結合論における $\sigma$ 結合と $\pi$ 結合の数理的な起源が明確に読み取れます（原子が z軸上に配置されている場合）。

$\sigma$ 結合 ($p_z - p_z$):項: $(P_z - A_z)(P_z - B_z) + \frac{1}{2p}$

解説: 第1項（幾何因子）は、軌道が「向かい合っている」効果を表します。通常、この項は負の値（ローブの符号が逆向きの領域での重なり）を持ちますが、第2項の拡散項 $\frac{1}{2p}$ と競合します。近距離では大きな重なりを持ちます。

符号: 通常の定義（正のローブ同士が向き合う）では、座標系の取り方により計算値は負になることがありますが（$z$正方向と$z$負方向）、物理的には結合性です。

$\pi$ 結合 ($p_x - p_x$ および $p_y - p_y$):項: $0 \cdot 0 + \frac{1}{2p} = \frac{1}{2p}$

解説: 結合軸方向（z）への変位がないため、幾何学的シフト項はゼロになります。純粋にガウス関数の「横方向の広がり」を表す $\frac{1}{2p}$ 項のみが生き残り、これが $\pi$ 重なりの正体です。$\sigma$ 結合に比べて距離減衰が急激であることも、この式の性質から説明できます。

直交成分 ($p_x - p_y$ など):項: $0 \cdot 0 + 0 = 0$

解説: クロネッカーのデルタ $\delta_{xy}=0$ かつ、幾何項もゼロであるため、重なりは厳密にゼロになります。