Full_CI.py

import numpy as np
from numpy import linalg as LA
from itertools import combinations
import pyscf
from pyscf import gto, scf, fci

np.set_printoptions(precision=8, suppress=True)

def hydrogen_chain(n_atoms, distance=1.0):
    atom_str = ''
    for i in range(n_atoms):
        atom_str += f'H  {i * distance:.6f}  0.000000  0.000000\n'
    return atom_str

mol = gto.M(
    atom = hydrogen_chain(6, distance=1.0),
    basis = 'sto-3g', unit = 'Bohr'
)

mf = scf.RHF(mol); mf.kernel()
enuc = mol.energy_nuc()
E_HF = mf.e_tot

C = mf.mo_coeff
N_occ = mol.nelectron // 2
N_spatial = mol.nao
N_spin = 2 * N_spatial
N_elec = mol.nelectron

h_ao = mol.intor('int1e_kin') + mol.intor('int1e_nuc')
h_mo = C.T @ h_ao @ C
eri_ao = mol.intor('int2e')
eri_mo = np.einsum('uP,vQ,uvls,lR,sS->PQRS', C, C, eri_ao, C, C, optimize=True)

h_spin = np.zeros((N_spin, N_spin))
for p in range(N_spin):
    for q in range(N_spin):
        if (p % 2) == (q % 2):
            h_spin[p,q] = h_mo[p//2, q//2]

eri_phys = np.zeros((N_spin, N_spin, N_spin, N_spin))
for p in range(N_spin):
    for q in range(N_spin):
        for r in range(N_spin):
            for s in range(N_spin):
                if (p % 2) == (r % 2) and (q % 2) == (s % 2):
                    eri_phys[p,q,r,s] = eri_mo[p//2, r//2, q//2, s//2]

as_eri = eri_phys - eri_phys.transpose(0,1,3,2)

dets = list(combinations(range(N_spin), N_elec))
n_dets = len(dets)
print(f"Number of determinants: {n_dets}")

# Slater-Condon rules
def sign_and_alignment(d1, d2):
    """Find the alignment (common orbitals, differences) and sign.
    Returns (common, unique1, unique2, sign) where common is sorted common 
    occupied, unique1 are in d1 not d2, unique2 in d2 not d1."""
    s1 = set(d1); s2 = set(d2)
    u1 = sorted(s1 - s2); u2 = sorted(s2 - s1)
    common = sorted(s1 & s2)
    # Compute sign from moving orbitals to align
    # Place unique of d1 at front, sorted; unique of d2 at front, sorted
    # Count swaps needed
    # For d1: move u1 orbitals to front (in order), count transpositions
    # d1 is ordered. Position of u1[k] in d1 is some index.
    def count_swaps(det, uniques):
        det = list(det)
        swaps = 0
        for k, u in enumerate(uniques):
            pos = det.index(u)
            # move to position k
            swaps += (pos - k)
            det.pop(pos); det.insert(k, u)
        return swaps
    n1 = count_swaps(d1, u1)
    n2 = count_swaps(d2, u2)
    sign = (-1)**(n1 + n2)
    return common, u1, u2, sign

def H_element(d1, d2):
    common, u1, u2, sign = sign_and_alignment(d1, d2)
    ndiff = len(u1)
    if ndiff == 0:
        # Same determinant
        h_sum = sum(h_spin[p,p] for p in d1)
        v_sum = 0.0
        for idx_p, p in enumerate(d1):
            for q in d1[idx_p+1:]:
                v_sum += as_eri[p,q,p,q]
        return h_sum + v_sum
    elif ndiff == 1:
        m = u1[0]; p = u2[0]
        val = h_spin[m,p]
        for n in common:
            val += as_eri[m,n,p,n]
        return sign * val
    elif ndiff == 2:
        m,n = u1; p,q = u2
        return sign * as_eri[m,n,p,q]
    else:
        return 0.0

print("Building FCI Hamiltonian matrix...")
H_FCI = np.zeros((n_dets, n_dets))
for I in range(n_dets):
    for J in range(I, n_dets):
        val = H_element(dets[I], dets[J])
        H_FCI[I,J] = val
        H_FCI[J,I] = val

print("Diagonalizing...")
eigvals, eigvecs = np.linalg.eigh(H_FCI)

print(f"\nHF total energy:   {E_HF:.8f} Ha")
print(f"FCI ground state:  {eigvals[0] + enuc:.8f} Ha")
print(f"FCI correlation:   {eigvals[0] + enuc - E_HF:.8f} Ha")

cisolver = fci.FCI(mf)
E_fci_pyscf, ci = cisolver.kernel()
print(f"\nPySCF FCI:         {E_fci_pyscf:.8f} Ha")
print(f"PySCF correlation: {E_fci_pyscf - E_HF:.8f} Ha")
print(f"Match: {abs(eigvals[0] + enuc - E_fci_pyscf) < 1e-8}")

print(f"\nLowest 5 FCI states (correlation energy):")
for n in range(5):
    print(f"  State {n}: E = {eigvals[n] + enuc:.6f} Ha, E_corr = {eigvals[n] + enuc - E_HF:.6f} Ha")

print("\n--- Full CI with Davidson-Liu ---")
print(f"Number of determinants: {n_dets}")

H_diag = np.array([H_element(dets[I], dets[I]) for I in range(n_dets)])

def sigma(v):
    s = np.zeros(n_dets)
    for I in range(n_dets):
        for J in range(n_dets):
            hij = H_element(dets[I], dets[J])
            s[I] += hij * v[J]
    return s

max_d_l_iter = 50
max_d_l_subspace = 20
tol = 1e-8
n_roots = 1

hf_det = tuple(range(N_elec))
hf_idx = dets.index(hf_det)
b = np.zeros(n_dets)
b[hf_idx] = 1.0

B = [b]
S = []

print("Starting Davidson-Liu iterations...")
for iteration in range(max_d_l_iter):
    sig = sigma(B[-1])
    S.append(sig)

    n_sub = len(B)
    H_sub = np.zeros((n_sub, n_sub))
    for i in range(n_sub):
        for j in range(n_sub):
            H_sub[i, j] = B[i] @ S[j]

    eigvals_sub, eigvecs_sub = LA.eigh(H_sub)

    E_d_l = eigvals_sub[0]
    c = eigvecs_sub[:, 0]

    x = sum(c[k] * B[k] for k in range(n_sub))
    Hx = sum(c[k] * S[k] for k in range(n_sub))
    r = Hx - E_d_l * x

    r_norm = LA.norm(r)
    print(f" Iter {iteration+1:2d}: E = {E_d_l:.10f}, |r| = {r_norm:.2e}, subspace = {n_sub}")

    if r_norm < tol:
        print(f"\nDavidson-Liu converged at iteration {iteration+1}")
        break

    t = np.zeros(n_dets)
    for I in range(n_dets):
        denom = E_d_l - H_diag[I]
        if abs(denom) > 1e-12:
            t[I] = -r[I] / denom
        else:
            t[I] = -r[I]

    for bk in B:
        t -= np.dot(t, bk) * bk
    norm_t = np.linalg.norm(t)
    if norm_t < 1e-14:
        print("  New vector too small, stopping")
        break
    t /= norm_t
    B.append(t)

    if len(B) > max_d_l_subspace:
        print(f"  Restarting subspace (was {len(B)})")
        x_restart = sum(c[k] * B[k] for k in range(n_sub))
        x_restart /= np.linalg.norm(x_restart)
        B = [x_restart]
        S = [sigma(x_restart)]


E_FCI_dav = E_d_l + enuc
print(f"\nFCI energy (Davidson): {E_FCI_dav:.8f} Ha")
print(f"FCI correlation:       {E_FCI_dav - (E_HF):.8f} Ha")