last_modified: 2026-01-23

注記: d軌道には、5成分の「球面調和関数型（Spherical Harmonics, $5d$）」と、6成分の「デカルト型（Cartesian, $6d$）」が存在する。計算機実装においては、積分の容易さから**Cartesian GTO（6d）**を計算し、必要に応じて変換行列を用いてSpherical（5d）へ変換するのが標準的である。本稿ではCartesian型（$d_{xx}, d_{xy}, \dots$）を扱う。 注記2: 生成AIへの指示で自動生成し、検算をしておりますが、内容が誤っている可能性があります。ご了承ください。Gemini 3.0 Pro及びGPT-5.2の性能を確かめる目的も兼ねています。

1. 数理的定式化#

中心 $A$ の短縮s軌道 $\Psi_A$ と、中心 $B$ の短縮d軌道成分 $\Psi_{B, \gamma\delta}$（$\gamma, \delta \in \{x, y, z\}$）の重なり積分 $S_{\gamma\delta}$ は以下の通り定義される。

$S_{\gamma\delta} = \sum_{\mu=1}^{K_A} \sum_{\nu=1}^{K_B} \underbrace{(d_{A,\mu} N_{A,\mu}^{prim})}_{D_{A,\mu}} \underbrace{(d_{B,\nu} N_{B,\nu}^{prim}(\gamma\delta))}_{D_{B,\nu}(\gamma\delta)} \cdot \langle \phi_{s,\mu} | \phi_{d_{\gamma\delta},\nu} \rangle_{un-norm}$

ここで、数式を「正規化係数部分」と「未正規化積分部分」に明確に分離する。

1.1 プリミティブ正規化定数 $N^{prim}$#

Cartesian GTO $\phi(\mathbf{r}) = x^l y^m z^n e^{-\zeta r^2}$

の正規化定数は、角運動量量子数 $\mathbf{L}=(l,m,n)$ に依存する。

s軌道 (0,0,0): $N_s = (2\zeta/\pi)^{3/4}$

d軌道 $xx$ (2,0,0): $N_{xx} = (2\zeta/\pi)^{3/4} \sqrt{ (8\zeta)^2 / 3 } = \frac{8\zeta}{\sqrt{3}} (2\zeta/\pi)^{3/4}$

d軌道 $xy$ (1,1,0): $N_{xy} = (2\zeta/\pi)^{3/4} \sqrt{ (8\zeta)^2 / 1 } = 8\zeta (2\zeta/\pi)^{3/4}$

重要: $N_{xy} = \sqrt{3} N_{xx}$ の関係があり、ここを区別しないと積分値が $\sqrt{3}$ 倍ずれる。

1.2 未正規化重なり積分 $\langle \phi | \phi \rangle_{un-norm}$#

正規化定数を除いた、純粋なガウス積分部分は以下の漸化式で表される（Obara-Saikaの初項）。ここで $S_{ss}^{(0)}$ は未正規化s-s重なり積分とする。

これを用いると、s-d重なり積分の「積分部分」は次の一元化公式となる。

$\langle s | d_{\gamma\delta} \rangle_{un-norm} = S_{ss}^{(0)} \cdot \left[ \mathcal{G}_\gamma \mathcal{G}_\delta + \frac{\delta_{\gamma\delta}}{2(\alpha_\mu + \beta_\nu)} \right]$

$\mathcal{G}_\gamma = \frac{\alpha_\mu (A_\gamma - B_\gamma)}{\alpha_\mu + \beta_\nu}$

$\delta_{\gamma\delta}$: クロネッカーのデルタ（$\gamma=\delta$ の時 $1$、それ以外 $0$）

2. Pythonによる実装#

import numpy as np

# --- 1. 数学ユーティリティ ---

def double_factorial(n):
    """二重階乗 n!! (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:
    """
    Cartesian GTOプリミティブの正規化定数を計算する (厳密解)。
    N = (2zeta/pi)^0.75 * sqrt( (8zeta)^L / (prod (2l-1)!!) )
    """
    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(alpha: float, beta: float, 
                                             lmn_A: tuple, lmn_B: tuple, 
                                             R_vec: np.ndarray) -> float:
    """
    一般化された未正規化プリミティブ重なり積分 (今回はs-d用およびd-d正規化用)。
    R_vec = A - B
    ※ここでは本題の s-d 公式と、d-d 正規化に必要な R=0 での積分を兼用する簡易ロジック
    """
    p = alpha + beta
    R2 = np.dot(R_vec, R_vec)
    
    # 基本のs-s積分 (未正規化)
    sss_0 = (np.pi / p)**1.5 * np.exp( - (alpha * beta / p) * R2 )
    
    # Aがs軌道 (0,0,0) かつ Bがd軌道 (lx,ly,lz) の場合
    if lmn_A == (0,0,0) and sum(lmn_B) == 2:
        lx, ly, lz = lmn_B
        # Gベクトル: G = alpha * (A - B) / p
        G = (alpha / p) * R_vec 
        
        # d_xx (2,0,0) -> Gx^2 + 1/2p
        if lx == 2: return sss_0 * (G[0]**2 + 1.0/(2.0*p))
        if ly == 2: return sss_0 * (G[1]**2 + 1.0/(2.0*p))
        if lz == 2: return sss_0 * (G[2]**2 + 1.0/(2.0*p))
        # d_xy (1,1,0) -> Gx*Gy
        if lx==1 and ly==1: return sss_0 * (G[0]*G[1])
        if lx==1 and lz==1: return sss_0 * (G[0]*G[2])
        if ly==1 and lz==1: return sss_0 * (G[1]*G[2])

    # 正規化計算用: 中心が同じ(R=0) で d-d 重なり (同じ成分同士)
    # Integral( x^2N * exp(-p r^2) ) の形
    if np.allclose(R2, 0.0) and lmn_A == lmn_B:
        # 例: <xx|xx> -> x^4項。 <xy|xy> -> x^2 y^2項。
        # 一般公式: I = prod( (2(li+lj)-1)!! / (2p)^(li+lj) ) * (pi/p)^1.5
        # ここでは lmn_A = lmn_B = (lx, ly, lz) なので power terms は 2lx, 2ly, 2lz
        lx, ly, lz = lmn_A
        
        # x成分の積分: (2*(2lx)-1)!! / (4p)^lx * sqrt(pi/p) ... いやもっと単純に
        # int x^2k exp(-p x^2) = (2k-1)!! / (2p)^k * sqrt(pi/p)
        # 今回 k = 2lx (合成した多項式の次数) ではなく、xの次数合計は 2*lx
        
        term_x = double_factorial(2*(2*lx) - 1) / ((4.0*p)**lx) # ? 少し複雑
        # 正規化用なので解析解をハードコードで安全にいく
        # <xx|xx> (total x^4) -> 3!! / (2p)^2 * sqrt(pi/p) * sqrt(pi/p)^2
        # <xy|xy> (total x^2 y^2) -> 1!!/(2p) * 1!!/(2p) * ...
        
        # 汎用1D積分 int t^n exp(-p t^2) dt = (n-1)!! / (2p)^(n/2) * sqrt(pi/p) (n:even)
        def int_1d(n, p_val):
            return double_factorial(n-1) / ((2.0*p_val)**(n/2)) * np.sqrt(np.pi/p_val)
            
        Ix = int_1d(2*lx, p) # xの次数は lx + lx = 2lx
        Iy = int_1d(2*ly, p)
        Iz = int_1d(2*lz, p)
        return Ix * Iy * Iz

    return 0.0 # 今回のスコープ外

# --- 2. 短縮係数の正規化処理 (重要) ---

def get_contraction_normalizer(exponents: list, coeffs: list, lmn: tuple) -> float:
    """
    短縮GTOの正規化定数 N_cont を計算する。
    1 = N_cont^2 * sum_ij [ c_i c_j <phi_i | phi_j> ]
    """
    n_prim = len(exponents)
    overlap_sum = 0.0
    
    for i in range(n_prim):
        for j in range(n_prim):
            # 1. プリミティブ正規化定数の積
            N_prim_i = get_primitive_norm(exponents[i], *lmn)
            N_prim_j = get_primitive_norm(exponents[j], *lmn)
            
            # 2. 未正規化プリミティブ間の重なり (R=0, 同一成分)
            # alpha_i, alpha_j が異なる場合もここで厳密に計算される
            S_prim_unnorm = calculate_unnormalized_primitive_overlap(
                exponents[i], exponents[j], lmn, lmn, np.zeros(3)
            )
            
            overlap_sum += coeffs[i] * coeffs[j] * N_prim_i * N_prim_j * S_prim_unnorm
            
    return 1.0 / np.sqrt(overlap_sum)

# --- 3. s-d 重なり積分計算クラス (修正版) ---

def calculate_overlap_sd_rigorous(
    center_A: np.ndarray, exps_A: list, coeffs_A: list, # s-orbital
    center_B: np.ndarray, exps_B: list, coeffs_B: list  # d-orbital
) -> dict:
    
    R_vec = center_A - center_B
    d_components = {
        'xx': (2, 0, 0), 'yy': (0, 2, 0), 'zz': (0, 0, 2),
        'xy': (1, 1, 0), 'xz': (1, 0, 1), 'yz': (0, 1, 1)
    }
    
    # 1. A(s) の正規化
    N_cont_A = get_contraction_normalizer(exps_A, coeffs_A, (0,0,0))
    
    # 2. B(d) の正規化 (成分ごとに計算！)
    # 通常の量子化学計算ではd殻で1つの係数を使うが、Cartesianでは成分でNormが変わる
    N_cont_B_dict = {}
    for label, lmn in d_components.items():
        N_cont_B_dict[label] = get_contraction_normalizer(exps_B, coeffs_B, lmn)

    results = {k: 0.0 for k in d_components}

    # 3. 二重ループ計算
    for i, alpha in enumerate(exps_A):
        # 係数 × プリミティブ正規化定数 (A)
        # N_prim は指数ごとに異なる
        D_A = coeffs_A[i] * N_cont_A * get_primitive_norm(alpha, 0,0,0)
        
        for j, beta in enumerate(exps_B):
            
            # 共通部分の計算
            p = alpha + beta
            sss_0 = (np.pi / p)**1.5 * np.exp( - (alpha * beta / p) * np.dot(R_vec, R_vec) )
            G = (alpha / p) * R_vec
            inv_2p = 1.0 / (2.0 * p)
            
            for label, (lx, ly, lz) in d_components.items():
                # 係数 × プリミティブ正規化定数 (B)
                # ここで beta と (lx,ly,lz) に依存する N_prim を正確に掛ける
                D_B = coeffs_B[j] * N_cont_B_dict[label] * get_primitive_norm(beta, lx, ly, lz)
                
                # 重なり項 (Geometric Factor)
                term = 0.0
                if lx == 2:   term = G[0]**2 + inv_2p  # xx
                elif ly == 2: term = G[1]**2 + inv_2p  # yy
                elif lz == 2: term = G[2]**2 + inv_2p  # zz
                elif lx==1 and ly==1: term = G[0]*G[1] # xy
                elif lx==1 and lz==1: term = G[0]*G[2] # xz
                elif ly==1 and lz==1: term = G[1]*G[2] # yz
                
                # 総和: D_A * D_B * S_unnorm
                results[label] += D_A * D_B * sss_0 * term

    return results

# --- 検証 ---
if __name__ == "__main__":
    # テストケース: x軸上に配置
    pos_A = np.array([0.0, 0.0, 0.0])
    pos_B = np.array([2.0, 0.0, 0.0])
    
    # 簡易な係数 (正規化定数の動作確認のため、同じ指数・係数を入れる)
    exps = [1.0]
    coeffs = [1.0]
    
    overlaps = calculate_overlap_sd_rigorous(pos_A, exps, coeffs, pos_B, exps, coeffs)
    
    print("Rigorous s-d Overlap Results:")
    for k, v in overlaps.items():
        print(f"  {k}: {v:.6f}")
        
    print("\nCheck Normalization Consistency:")
    # xx と xy の正規化が正しく機能していれば、
    # 物理的に等価な配置（回転して重なるようなケース）で整合性が取れるはずである。
    # ここでは、xxが大きく、yy/zzが小さく(inv_2p項のみ)、混合項が0になることを確認。

3. 物理的解釈の注意点：6d vs 5d#

計算結果の解釈において、以下の点に注意が必要である。

3.1 Cartesian (6d) 特有の挙動#

上記の実装は「6成分 Cartesian d軌道」に基づいている。$\{ x^2, y^2, z^2, xy, xz, yz \}$この基底では、s軌道（球対称）との重なりにおいて、$d_{yy}$ や $d_{zz}$ 成分もゼロにならない（値 $\propto \frac{1}{2p} S_{ss}$）。これは Cartesian d軌道の $y^2, z^2$ 成分が、球対称成分（$r^2$）を含んでいるためである（s軌道コンタミネーション）。

3.2 Spherical (5d) への変換#

物理的に純粋なd軌道（$l=2$ の固有状態）を議論する場合、以下の5成分への変換が必要となる。特に、$d_{z^2}$（あるいは $d_{3z^2-r^2}$）は以下のように構成される。$d_{z^2}^{(5d)} \propto 2 d_{zz}^{(6d)} - d_{xx}^{(6d)} - d_{yy}^{(6d)}$もし $s$ 軌道と $d$ 軌道が $z$ 軸上で対向している場合、Spherical 表記では $d_{z^2}$ のみが $s$ 軌道と重なりを持ち、$\pi$ 対称性を持つ $d_{xz}, d_{yz}$ や、$\delta$ 対称性を持つ $d_{xy}, d_{x^2-y^2}$ は厳密にゼロとなる。Cartesian 計算結果を解析に使用する際は、この変換を意識する必要がある。