CI計算における多重項の縮退数の不一致

import numpy as np
import pandas as pd
import itertools
from math import comb

# ========= ユーザ設定 =========
N_ORB = 10     # 基底関数数
N_ELEC = 5     # 電子数
ONE_PKL = 'one_integral.pkl'  # columns: ['i','j','one_integral']（1始まり）
TWO_PKL = 'two_integral.pkl'  # columns: ['i','j','k','l','integral']（1始まり）

# =============================

# ---- 基底（スレーター行列式＝占有0/1ベクトル）を生成 ----
def generate_slater_basis(n_orb, n_elec):
    basis_bits = []
    basis_occ_lists = []
    for occ in itertools.combinations(range(n_orb), n_elec):
        bits = np.zeros(n_orb, dtype=np.int8)
        bits[list(occ)] = 1
        basis_bits.append(bits)
        basis_occ_lists.append(list(occ))  # 0始まりの占有位置
    return basis_bits, basis_occ_lists

B_BITS, B_OCC = generate_slater_basis(N_ORB, N_ELEC)   #長さ=252

# ---- 反交換に基づく位相：左から s-1 番目までの占有数の偶奇 ----
def parity_before(r, s):
    # r: 0/1 の occupation（0始まり）
    # s: 0始まりの軌道インデックス
    return (-1) ** int(r[:s].sum())

# ---- ⟨p| a_i† a_j |q⟩（i,j は 0始まり）----
def me_1body(p_bits, q_bits, i, j):
    # a_j が q に作用できる条件
    if q_bits[j] == 0:
        return 0
    # a_i† が最終的に作用できる条件
    if p_bits[i] == 0:
        return 0

    # 逐次作用して p に一致するか確認
    r = q_bits.copy()
    # a_j
    phase_j = parity_before(r, j)
    r[j] = 0
    # a_i†
    if r[i] == 1:  # Pauli
        return 0
    phase_i = parity_before(r, i)
    r[i] = 1

    if np.array_equal(r, p_bits):
        if i == j:
            # a_i† a_i = n_i → 対角で 1、位相は (+1)*(+1) になる
            # ただし上の条件で p_bits[i]==1 かつ q_bits[i]==1 なのでOK
            return 1
        return phase_j * phase_i
    return 0

# ---- ⟨p| a_i† a_j† a_l a_k |q⟩（i,j,k,l は 0始まり）----
def me_2body(p_bits, q_bits, i, j, k, l):
    # 同一演算子の二連はゼロ
    if i == j or k == l:
        return 0
    # 消滅は q の占有から
    if q_bits[k] == 0 or q_bits[l] == 0:
        return 0
    # 生成は最終 p の占有へ
    if p_bits[i] == 0 or p_bits[j] == 0:
        return 0

    # 逐次作用：a_k → a_l → a_j† → a_i†（右から左）
    r = q_bits.copy()

    # a_k
    phase1 = parity_before(r, k)
    r[k] = 0

    # a_l
    phase2 = parity_before(r, l)
    if r[l] == 0:
        return 0
    r[l] = 0

    # a_j†
    phase3 = parity_before(r, j)
    if r[j] == 1:
        return 0
    r[j] = 1

    # a_i†
    phase4 = parity_before(r, i)
    if r[i] == 1:
        return 0
    r[i] = 1

    if np.array_equal(r, p_bits):
        return phase1 * phase2 * phase3 * phase4
    return 0

# ---- 積分を読み込み（1始まり） → ルックアップ辞書に ----
one_df = pd.read_pickle(ONE_PKL)
two_df = pd.read_pickle(TWO_PKL)
print(one_df.shape)
print(two_df.shape)


one_lookup = {
    (int(i), int(j)): complex(np.real(v)) if pd.notna(v) else 0.0
    for i, j, v in one_df[['i', 'j', 'one_integral']].itertuples(index=False)
}

two_lookup = {
    (int(i), int(j), int(k), int(l)): complex(np.real(v)) if pd.notna(v) else 0.0
    for i, j, k, l, v in two_df[['i', 'j', 'k', 'l', 'integral']].itertuples(index=False)
}

# ---- CI 行列の構築（252×252）----
NB = len(B_BITS)
H = np.zeros((NB, NB), dtype=complex)

# 励起次数で枝刈り：1体は最大1電子移動、2体は最大2電子移動
def hamming(x, y):
    return int(np.sum(x ^ y))  # xor の1の個数（np.int8でもOK）

for p in range(NB):
    p_bits = B_BITS[p]
    p_occ = B_OCC[p]  # 0始まり
    # p側で 1体：i は pで占有の軌道（0始まり）
    p_occ_1based = [i+1 for i in p_occ]
    # p側で 2体：i<j の2個選び（0始まり）
    p_pairs_0 = list(itertools.combinations(p_occ, 2))
    p_pairs_1b = [(i+1, j+1) for (i,j) in p_pairs_0]

    for q in range(NB):
        q_bits = B_BITS[q]
        q_occ = B_OCC[q]
        q_occ_1based = [k+1 for k in q_occ]
        # 励起次数（xorの1の数）は 0,2,4 のみ寄与
        d = hamming(p_bits, q_bits)
        if d not in (0, 2, 4):
            continue

        # ---- 一電子項：∑_{ij} h_ij ⟨p|a_i† a_j|q⟩ ----
        # 条件：j は q で占有、i は p で占有（i,j は 0始まり）
        term1 = 0.0
        for i0 in p_occ:        # 0始まり
            for j0 in q_occ:    # 0始まり
                val = me_1body(p_bits, q_bits, i0, j0)
                if val == 0:
                    continue
                hij = one_lookup.get((i0+1, j0+1), 0.0)  # 1始まりで検索
                if hij != 0.0:
                    term1 += val * hij

        # ---- 二電子項： 1/2 ∑_{ijkl} (ij|kl) ⟨p|a_i† a_j† a_l a_k|q⟩ ----
        term2_sum = 0.0
        # i<j, k<l で十分だが、位相は me_2body に任せるので全順序にしてもOK。
        # ただし計算量を抑えるため i<j, k<l で回す。
        for (i0, j0) in p_pairs_0:   # 0始まり
            for (k0, l0) in itertools.combinations(q_occ, 2):
                # 4通りの順序を試す（(i,j)の順、(k,l)の順）※辞書は1始まり
                # Coulomb integral は (ij|kl) = (ji|lk) = (kl|ij) ... 対称だが、
                # two_lookup が完全表でない可能性に備えて必要な順序を全部問い合わせる。
                contrib = 0.0
                # (i,j,k,l)
                for (ii,jj) in ((i0,j0),(j0,i0)):
                    for (kk,ll) in ((k0,l0),(l0,k0)):
                        val = me_2body(p_bits, q_bits, ii, jj, kk, ll)
                        if val == 0:
                            continue
                        gij = two_lookup.get((ii+1, jj+1, kk+1, ll+1), 0.0)
                        if gij != 0.0:
                            contrib += val * gij
                term2_sum += contrib

        H[p, q] = term1 + 0.5 * term2_sum

# 数値非対称を消す（理論上 H は対称）
H = 0.5 * (H + H.T)

# ---- 対角化 ----
eigvals, eigvecs = np.linalg.eigh(H)

print("Matrix size:", H.shape)
print(eigvals[:40])