Crawford7_pyscf.py

import numpy as np
from numpy import linalg as LA
import time
import pyscf
from pyscf import gto
from pyscf import lo
from itertools import combinations

np.set_printoptions(precision=7, suppress=True, linewidth=150)

mol = gto.M(
    atom = '''
    O   0.000000000000  -0.143225816552   0.000000000000
    H   1.638036840407   1.136548822547  -0.000000000000
    H  -1.638036840407   1.136548822547  -0.000000000000
    ''',
    basis = 'sto-3g',
    unit = 'Bohr'
)

enuc = mol.energy_nuc()

dimension = mol.nao
N_elec = mol.nelectron
N_occ = mol.nelectron // 2
s_matrix, t_matrix, v_matrix = np.zeros((dimension, dimension), dtype=float), np.zeros((dimension, dimension), dtype=float), np.zeros((dimension, dimension), dtype=float)
eri_tensor = np.zeros((dimension, dimension, dimension, dimension))
delta_threshold = 1.0e-8
rms_threshold = 1.0e-5

s_matrix = mol.intor('int1e_ovlp')
t_matrix = mol.intor('int1e_kin')
v_matrix = mol.intor('int1e_nuc')

H_core = t_matrix + v_matrix

eri_tensor = mol.intor('int2e')

def diagonalization_S(S):
    eigvals,eigvecs = LA.eigh(S)
    Lambda_S = np.diag((eigvals) ** (-0.5))
    X = np.dot(np.dot(eigvecs, Lambda_S), eigvecs.T)
    return X

def density(X,F):
    F_prime = np.dot(np.dot(X.T, F), X)
    eigvals,eigvecs = LA.eigh(F_prime)
    C = np.dot(X, eigvecs)
    C_occ = C[:, :N_occ]
    D = np.dot(C_occ, C_occ.T)
    return D, C, eigvals

def scf(D, H, F):
    return np.sum(D * (H + F))

def fock(D, H, eri):
    F = H.copy()
    J = np.einsum('pqrs,rs->pq', eri, D)
    K = np.einsum('prqs,rs->pq', eri, D)
    F = H + 2*J - K
    return F

def error_matrix(F, D):
    term1 = np.dot(np.dot(F,D),s_matrix)
    term2 = np.dot(np.dot(s_matrix,D),F)
    return term1-term2

def diis_equation(error_matrices):
    m = len(error_matrices)
    B = np.zeros((m+1, m+1))

    for i in range(m):
        for j in range(m):
            B[i,j] = np.sum(error_matrices[i] * error_matrices[j])
    B[-1,:-1] = -1.0
    B[:-1,-1] = -1.0
    B_vector = np.zeros(m+1)
    B_vector[-1] = -1.0

    c_vector = np.linalg.solve(B,B_vector)
    return c_vector[:-1]

def diis_fock(c_coeff, fock_matrices):
    fock_stack = np.stack(fock_matrices, axis=0)
    return np.einsum('i, ijk -> jk', c_coeff, fock_stack, optimize=True)
    

MAX_DIIS = 6
diis_fock_matrices = []
diis_error_matrices = []
X = diagonalization_S(s_matrix)
D, C, orbital_energies = density(X, H_core)
E_elec_old = None
D_old = D.copy()

time_hf_start = time.time()
# Iter 0: initial guess energy (Crawford convention)
E_elec = scf(D, H_core, H_core)
print(f"Circle 0, E_elec = {E_elec} Hartree, E_tol = {E_elec + enuc} Hartree.")
E_elec_old = E_elec
 
for i in range(1, 100):
    # 1. Build Fock from current D
    F = fock(D, H_core, eri_tensor)
 
    # 2. Energy with matched (D, F)
    E_elec = scf(D, H_core, F)
 
    # 3. Error with the same (F, D) pair
    error = error_matrix(F, D)
 
    # 4. Store matched (F, error) into DIIS
    diis_fock_matrices.append(F.copy())
    diis_error_matrices.append(error.copy())
    if len(diis_fock_matrices) > MAX_DIIS:
        diis_fock_matrices.pop(0)
        diis_error_matrices.pop(0)
 
    # 5. Convergence check
    delta_E = np.abs(E_elec - E_elec_old)
    rms_D = np.sqrt(np.sum((D - D_old)**2))
    print(f"Circle {i}, E_elec = {E_elec} Hartree, E_tol = {E_elec + enuc} Hartree.")
    print(f"Convergence Check: Delta_E = {delta_E:.2e}, RMS_D = {rms_D:.2e}")
    if delta_E < delta_threshold and rms_D < rms_threshold:
        print("\n--- SCF Converged ---")
        print("Final Total Energy (E_tol):", E_elec + enuc)
 
        print("\nOrbital Energies (eigvals from diagonalization):")
        print(orbital_energies)
 
        print("\nPerforming AO-to-MO transformation: F_mo = C.T @ F @ C")
 
        F_mo = C.T @ F @ C
 
        print("\nMO-Basis Fock Matrix (F_mo):")
        print(F_mo)
 
        print("\nMO coefficients matrix:")
        print(C)
        break
 
    # 6. DIIS extrapolation
    if len(diis_fock_matrices) > 1:
        F_extrapolated = diis_fock(diis_equation(diis_error_matrices), diis_fock_matrices)
    else:
        F_extrapolated = F
 
    # 7. New D from extrapolated Fock
    D_old = D.copy()
    D, C, orbital_energies = density(X, F_extrapolated)
    E_elec_old = E_elec


print("\n--- Starting MP2 Calculation ---")
print("Preparing the MO basis ERI...")

step_1 = np.einsum('uvls, sS -> uvlS', eri_tensor, C, optimize=True)
step_2 = np.einsum('uvlS, lR -> uvRS', step_1, C, optimize=True)
step_3 = np.einsum('uvRS, vQ -> uQRS', step_2, C, optimize=True)
eri_tensor_MO = np.einsum('uQRS, uP -> PQRS', step_3, C, optimize=True)
print(f"MO ERI tensor (pq|rs) created with shape: {eri_tensor_MO.shape}")

print("\nPreparing the MP2 energy...")
eps_occ = orbital_energies[:N_occ]
eps_vir = orbital_energies[N_occ:]

eps_i = eps_occ[:,None,None,None]
eps_j = eps_occ[None,None,:,None]
eps_a = eps_vir[None,:,None,None]
eps_b = eps_vir[None,None,None,:]

denominator = eps_i + eps_j - eps_a - eps_b

O = slice(None,N_occ)
V = slice(N_occ,None)

I_iajb = eri_tensor_MO[O,V,O,V]
I_ibja = eri_tensor_MO[O,V,O,V].transpose(0, 3, 2, 1)

numerator = I_iajb * (2*I_iajb-I_ibja)

E_MP2_tensor = numerator / denominator
E_MP2 = np.sum(E_MP2_tensor)
print(f"\n--- MP2 Calculation Complete ---")
print(f"MP2 Correlation Energy: {E_MP2:.8f} Hartree")
print(f"Total HF Energy (E_tol):  {E_elec + enuc:.8f} Hartree")
print(f"Total MP2 Energy:         {E_elec + enuc + E_MP2:.8f} Hartree")

print("\n--- Starting CCSD Calculation ---")
print("Preparing the spin-orbital basis Fock Matrix...")
N_spatial = dimension
N_spin = N_spatial * 2
f_pq = np.diag(np.repeat(orbital_energies,2))
spin_orbital_energies = np.diag(f_pq)

print("Preparing spin-orbital eri tensor...")
eri_tensor_spin = 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):
                P = p // 2
                Q = q // 2
                R = r // 2
                S = s // 2

                sigma_p = p % 2
                sigma_q = q % 2
                sigma_r = r % 2
                sigma_s = s % 2
                
                if (sigma_p == sigma_r) and (sigma_q == sigma_s):
                    eri_tensor_spin[p,q,r,s] = eri_tensor_MO[P,R,Q,S]

print("Preparing the initial guess of T1 and T2...")
N_occ_spin = N_occ * 2
N_vir_spin = N_spin - N_occ_spin

T1 = np.zeros((N_occ_spin, N_vir_spin))

O_spin = slice(None, N_occ_spin)
V_spin = slice(N_occ_spin, None)

eps_occ_spin = spin_orbital_energies[O_spin]
eps_vir_spin = spin_orbital_energies[V_spin]

eps_i_spin = eps_occ_spin[:,None,None,None]
eps_j_spin = eps_occ_spin[None,:,None,None]
eps_a_spin = eps_vir_spin[None,None,:,None]
eps_b_spin = eps_vir_spin[None,None,None,:]

D_ijab = eps_i_spin + eps_j_spin - eps_a_spin - eps_b_spin

I_ijab = eri_tensor_spin[O_spin,O_spin,V_spin,V_spin]
I_ijba = eri_tensor_spin[O_spin,O_spin,V_spin,V_spin].transpose(0,1,3,2)
as_ijab = I_ijab - I_ijba

T2 = as_ijab / D_ijab

def tau_intermediates(T1, T2):
    t1_t1 = np.einsum('ia, jb -> ijab', T1, T1, optimize=True)
    T1_product = t1_t1 - t1_t1.transpose(0,1,3,2)
    tau = T2 + T1_product
    tau_tilde = T2 + 0.5 * T1_product
    return tau, tau_tilde

eri_VOVO = eri_tensor_spin[V_spin,O_spin,V_spin,O_spin]
eri_VOVV = eri_tensor_spin[V_spin,O_spin,V_spin,V_spin]
eri_OVVV = eri_tensor_spin[O_spin,V_spin,V_spin,V_spin]
eri_OOVV = eri_tensor_spin[O_spin,O_spin,V_spin,V_spin]
eri_OOOV = eri_tensor_spin[O_spin,O_spin,O_spin,V_spin]
eri_OOOO = eri_tensor_spin[O_spin,O_spin,O_spin,O_spin]
eri_VVVV = eri_tensor_spin[V_spin,V_spin,V_spin,V_spin]
eri_OVVO = eri_tensor_spin[O_spin,V_spin,V_spin,O_spin]
eri_OOVO = eri_tensor_spin[O_spin,O_spin,V_spin,O_spin]
eri_OVOV = eri_tensor_spin[O_spin,V_spin,O_spin,V_spin] 
eri_VVVO = eri_tensor_spin[V_spin,V_spin,V_spin,O_spin] 
eri_VVOV = eri_tensor_spin[V_spin,V_spin,O_spin,V_spin] 
eri_OVOO = eri_tensor_spin[O_spin,V_spin,O_spin,O_spin]

f_oo = f_pq[O_spin,O_spin]
f_vv = f_pq[V_spin,V_spin]
f_ov = f_pq[O_spin,V_spin]
D_ia = np.diag(f_oo)[:,None]-np.diag(f_vv)[None,:]

def F_intermediates(T1, tau_tilde):
    term_ae_1 = np.einsum('mf, mafe -> ae', T1, eri_OVVV, optimize=True)
    term_ae_2 = np.einsum('mf, maef -> ae', T1, eri_OVVV, optimize=True)
    term_ae_a = term_ae_1 - term_ae_2

    term_ae_3 = np.einsum('mnaf, mnef -> ae', tau_tilde, eri_OOVV, optimize=True)
    term_ae_4 = np.einsum('mnaf, nmef -> ae', tau_tilde, eri_OOVV, optimize=True)
    term_ae_b = -0.5 * (term_ae_3 - term_ae_4)

    F_ae = term_ae_a + term_ae_b

    term_mi_1 = np.einsum('ne, mnie -> mi', T1, eri_OOOV, optimize=True)
    term_mi_2 = np.einsum('ne, nmie -> mi', T1, eri_OOOV, optimize=True)
    term_mi_a = term_mi_1 - term_mi_2

    term_mi_3 = np.einsum('inef, mnef -> mi', tau_tilde, eri_OOVV, optimize=True)
    term_mi_4 = np.einsum('inef, nmef -> mi', tau_tilde, eri_OOVV, optimize=True)
    term_mi_b = 0.5 * (term_mi_3 -term_mi_4)

    F_mi = term_mi_a + term_mi_b

    term_me_a = np.einsum('nf, mnef -> me', T1, eri_OOVV, optimize=True)
    term_me_b = np.einsum('nf, nmef -> me', T1, eri_OOVV, optimize=True)

    F_me = term_me_a - term_me_b

    return F_ae, F_mi, F_me

def W_intermediates(T1, T2, tau):
    term_mnij_a = eri_OOOO - eri_OOOO.transpose(1,0,2,3)

    term_mnij_1 = np.einsum('je, mnie -> mnij', T1, eri_OOOV, optimize=True)
    term_mnij_2 = np.einsum('je, nmie -> mnij', T1, eri_OOOV, optimize=True)
    term_mnij_3 = np.einsum('ie, mnje -> mnij', T1, eri_OOOV, optimize=True)
    term_mnij_4 = np.einsum('ie, nmje -> mnij', T1, eri_OOOV, optimize=True)
    term_mnij_b = term_mnij_1 - term_mnij_2 - (term_mnij_3 -term_mnij_4)

    term_mnij_5 = np.einsum('ijef, mnef -> mnij', tau, eri_OOVV, optimize=True)
    term_mnij_6 = np.einsum('ijef, nmef -> mnij', tau, eri_OOVV, optimize=True)
    term_mnij_c = (term_mnij_5 - term_mnij_6) / 4

    W_mnji = term_mnij_a + term_mnij_b + term_mnij_c

    term_abef_a = eri_VVVV - eri_VVVV.transpose(1,0,2,3)

    term_abef_1 = np.einsum('mb, amef -> abef', T1, eri_VOVV, optimize=True)
    term_abef_2 = np.einsum('mb, amfe -> abef', T1, eri_VOVV, optimize=True)
    term_abef_3 = np.einsum('ma, bmef -> abef', T1, eri_VOVV, optimize=True)
    term_abef_4 = np.einsum('ma, bmfe -> abef', T1, eri_VOVV, optimize=True)
    term_abef_b = term_abef_1 - term_abef_2 - (term_abef_3 -term_abef_4)

    term_abef_5 = np.einsum('mnab, mnef -> abef', tau, eri_OOVV, optimize=True)
    term_abef_6 = np.einsum('mnab, nmef -> abef', tau, eri_OOVV, optimize=True)
    term_abef_c = (term_abef_5 - term_abef_6) / 4

    W_abef = term_abef_a - term_abef_b + term_abef_c

    term_mbej_a = eri_OVVO - eri_OVOV.transpose(0,1,3,2)

    term_mbej_1 = np.einsum('jf, mbef -> mbej', T1, eri_OVVV, optimize=True)
    term_mbej_2 = np.einsum('jf, mbfe -> mbej', T1, eri_OVVV, optimize=True)
    term_mbej_b = term_mbej_1 -term_mbej_2

    term_mbej_3 = np.einsum('nb, mnej -> mbej', T1, eri_OOVO, optimize=True)
    term_mbej_4 = np.einsum('nb, nmej -> mbej', T1, eri_OOVO, optimize=True)
    term_mbej_c = term_mbej_3 - term_mbej_4

    as_OOVV = eri_OOVV - eri_OOVV.transpose(0,1,3,2)
    term_mbej_5 = 0.5 * np.einsum('jnfb, mnef -> mbej', T2, as_OOVV, optimize=True)
    term_mbej_6 = np.einsum('jf, nb, mnef -> mbej', T1, T1, as_OOVV, optimize=True)
    term_mbej_d = term_mbej_5 + term_mbej_6

    W_mbej = term_mbej_a + term_mbej_b - term_mbej_c - term_mbej_d

    return W_mnji, W_abef, W_mbej

def update_T1(T1, T2, F_ae, F_mi, F_me):
    RHS_T1 = f_ov.copy()
    RHS_T1 +=  np.einsum('ie, ae -> ia', T1, F_ae, optimize=True)
    RHS_T1 -=  np.einsum('ma, mi -> ia', T1, F_mi, optimize=True)
    RHS_T1 +=  np.einsum('imae, me -> ia', T2, F_me, optimize=True)

    RHS_T1 -= np.einsum('mf, maif -> ia', T1, eri_OVOV, optimize=True)

    RHS_T1 += np.einsum('mf, mafi -> ia', T1, eri_OVVO, optimize=True)
    
    as_OVVV = eri_OVVV - eri_OVVV.transpose(0,1,3,2) 
    RHS_T1 -= 0.5 * np.einsum('imef, maef -> ia', T2, as_OVVV, optimize=True)
    
    I_nmei = eri_OOVO 
    I_nmie = eri_OOOV.transpose(0,1,3,2) 
    as_OOOV_nmei = I_nmei - I_nmie
    
    RHS_T1 -= 0.5 * np.einsum('mnae, nmei -> ia', T2, as_OOOV_nmei, optimize=True)
    
    return RHS_T1

def update_T2(T1, T2, tau, F_ae, F_mi, F_me, W_mnji, W_abef, W_mbej):

    RHS_T2 = np.zeros_like(T2) 

    as_OOVV = eri_OOVV - eri_OOVV.transpose(0, 1, 3, 2)
    RHS_T2 += as_OOVV

    X_e = F_ae - 0.5 * np.einsum('mb, me -> be', T1, F_me, optimize=True)
    X_ijab = np.einsum('ijae, be -> ijab', T2, X_e, optimize=True) 
    RHS_T2 += X_ijab - X_ijab.transpose(0, 1, 3, 2) 

    Y_m = F_mi + 0.5 * np.einsum('je, me -> mj', T1, F_me, optimize=True)
    Y_ijab = np.einsum('imab, mj -> ijab', T2, Y_m, optimize=True)
    RHS_T2 -= (Y_ijab - Y_ijab.transpose(1, 0, 2, 3)) 
 
    RHS_T2 += 0.5 * np.einsum('mnab, mnij -> ijab', tau, W_mnji, optimize=True)

    RHS_T2 += 0.5 * np.einsum('ijef, abef -> ijab', tau, W_abef, optimize=True)

    Z_ijab = np.einsum('imae, mbej -> ijab', T2, W_mbej, optimize=True) 
    RHS_T2 += Z_ijab \
            - Z_ijab.transpose(1,0,2,3) \
            - Z_ijab.transpose(0,1,3,2) \
            + Z_ijab.transpose(1,0,3,2)

    as_OVVO = eri_OVVO - eri_OVOV.transpose(0,1,3,2)
    K_ijab = np.einsum('ie, ma, mbej -> ijab', T1, T1, as_OVVO, optimize=True)
    RHS_T2 -= (K_ijab - K_ijab.transpose(1,0,2,3) - K_ijab.transpose(0,1,3,2) + K_ijab.transpose(1,0,3,2)) 

    as_VVVO = eri_VVVO - eri_VVOV.transpose(0,1,3,2)
    L_ijab = np.einsum('ie, abej -> ijab', T1, as_VVVO, optimize=True)
    RHS_T2 += L_ijab - L_ijab.transpose(1,0,2,3) 

    as_OVOO = eri_OVOO - eri_OVOO.transpose(0, 1, 3, 2)
    M_ijab = np.einsum('ma, mbij -> ijab', T1, as_OVOO, optimize=True)
    RHS_T2 -= (M_ijab - M_ijab.transpose(0,1,3,2)) 

    return RHS_T2

def E_CCSD(T1, T2):
    as_OOVV = eri_OOVV - eri_OOVV.transpose(0,1,3,2)
    term_1 = np.einsum('ijab, ijab ->', T2, as_OOVV, optimize=True)
    term_2 = np.einsum('ia, jb, ijab ->', T1, T1, as_OOVV, optimize=True)
    return (term_1 / 4 + term_2 / 2)

print("\n--- Starting CCSD Iteration Loop ---")

max_iter = 50
E_threshold = 1.0e-7
E_ccsd_old = 0.0

for i in range(max_iter):

    T1_old = T1.copy()
    T2_old = T2.copy()

    tau, tau_tilde = tau_intermediates(T1,T2)
    F_ae, F_mi, F_me = F_intermediates(T1,tau_tilde)
    W_mnji, W_abef, W_mbej = W_intermediates(T1,T2,tau)

    RHS_T1 = update_T1(T1, T2, F_ae, F_mi, F_me)
    T1_new = RHS_T1 / D_ia

    RHS_T2 = update_T2(T1, T2, tau, F_ae, F_mi, F_me, W_mnji, W_abef, W_mbej)
    T2_new = RHS_T2 / D_ijab

    E_ccsd_new = E_CCSD(T1_new,T2_new)

    delta_E = np.abs(E_ccsd_new - E_ccsd_old)
    rms_T1 = np.sqrt(np.mean((T1_new - T1_old)**2))
    rms_T2 = np.sqrt(np.mean((T2_new - T2_old)**2))
    
    print(f"Iter: {i+1:2d}   E_corr: {E_ccsd_new:.12f}   Delta_E: {delta_E:.2e}   RMS_T1: {rms_T1:.2e}   RMS_T2: {rms_T2:.2e}")
    
    if delta_E < E_threshold:
        print("\n--- CCSD Converged! ---")
        break
    
    T1 = T1_new
    T2 = T2_new
    E_ccsd_old = E_ccsd_new


E_CCSD_total = E_elec + enuc + E_ccsd_new 
print(f"\nFinal CCSD Correlation Energy: {E_ccsd_new:.8f} Hartree")
print(f"Final Total CCSD Energy:     {E_CCSD_total:.8f} Hartree")

print("\n--- Starting (T) Calculation ---")
e_i = np.diag(f_oo)[:,None,None,None,None,None]
e_j = np.diag(f_oo)[None,:,None,None,None,None]
e_k = np.diag(f_oo)[None,None,:,None,None,None]
e_a = np.diag(f_vv)[None,None,None,:,None,None]
e_b = np.diag(f_vv)[None,None,None,None,:,None]
e_c = np.diag(f_vv)[None,None,None,None,None,:]

D_ijkabc = e_i + e_j + e_k - e_a - e_b - e_c

print("Building 'Disconnected' base tensor V_d...")
as_OOVV = eri_OOVV - eri_OOVV.transpose(0,1,3,2) 
V_d = np.einsum('ia, jkbc -> ijkabc', T1, as_OOVV, optimize=True)

print("Building 'Connected' base tensor V_c...")

as_VOVV = eri_VOVV - eri_VOVV.transpose(0,1,3,2) 
TermA = np.einsum('jkae, eibc -> ijkabc', T2, as_VOVV, optimize=True)


as_OVOO = eri_OVOO - eri_OVOO.transpose(0,1,3,2) 
TermB = -np.einsum('imbc, majk -> ijkabc', T2, as_OVOO, optimize=True)

V_c = TermA + TermB 

# P(p/qr)f(pqr) = f(pqr) - f(qpr) - f(rqp)
# D = P(i/jk) [ P(a/bc)[V] ]
# D = P(i/jk) [ V(abc) - V(bac) - V(cba) ]
# D = [V(i,jk,a,bc) - V(i,jk,b,ac) - V(i,jk,c,ba)]  (f(pqr))
#   - [V(j,ik,a,bc) - V(j,ik,b,ac) - V(j,ik,c,ba)]  (-f(qpr))
#   - [V(k,ji,a,bc) - V(k,ji,b,ac) - V(k,ji,c,ba)]  (-f(rqp))
print("Applying 9-term permutation P(i/jk)P(a/bc)...")

def apply_permutation(V):
    # V(i,j,k, a,b,c)
    V_abc = V
    # V(i,j,k, b,a,c)
    V_bac = V.transpose(0,1,2, 4,3,5)
    # V(i,j,k, c,b,a)
    V_cba = V.transpose(0,1,2, 5,4,3)
    
    # f(pqr) -> [V(abc) - V(bac) - V(cba)]
    term_i = V_abc - V_bac - V_cba

    # V(j,i,k, a,b,c)
    V_jik = V.transpose(1,0,2, 3,4,5)
    # V(j,i,k, b,a,c)
    V_jbk = V.transpose(1,0,2, 4,3,5)
    # V(j,i,k, c,b,a)
    V_jck = V.transpose(1,0,2, 5,4,3)

    # -f(qpr) -> -[V(jik,abc) - V(jik,bac) - V(jik,cba)]
    term_j = V_jik - V_jbk - V_jck

    # V(k,j,i, a,b,c)
    V_kji = V.transpose(2,1,0, 3,4,5)
    # V(k,j,i, b,a,c)
    V_kbi = V.transpose(2,1,0, 4,3,5)
    # V(k,j,i, c,b,a)
    V_kci = V.transpose(2,1,0, 5,4,3)

    # -f(rqp) -> -[V(kji,abc) - V(kji,bac) - V(kji,cba)]
    term_k = V_kji - V_kbi - V_kci
    
    # D = f(pqr) - f(qpr) - f(rqp)
    return term_i - term_j - term_k

D_d = apply_permutation(V_d)
D_c = apply_permutation(V_c)

print("Permutations complete.")

Numerator = D_c * (D_c + D_d)
E_T_tensor = Numerator / D_ijkabc

E_T = (1.0 / 36.0) * np.sum(E_T_tensor)

E_CCSD_T_total = E_CCSD_total + E_T

print(f"\n--- CCSD(T) Calculation Complete ---")
print(f"E(T) Correction:    {E_T:.12f} Hartree")
print(f"Final Total CCSD(T) Energy: {E_CCSD_T_total:.12f} Hartree")

print("\n--- Starting CIS Computation ---")
N_vir = dimension - N_occ

H_singlet = np.zeros((N_occ * N_vir, N_occ * N_vir))
H_triplet = np.zeros((N_occ * N_vir, N_occ * N_vir))

for i in range(N_occ):
    for a in range(N_vir):
        ia = i * N_vir + a
        for j in range(N_occ):
            for b in range(N_vir):
                jb = j * N_vir + b
                A = a + N_occ
                B = b + N_occ
                
                diag = (orbital_energies[A] - orbital_energies[i]) * (i==j) * (a==b)
                coulomb = eri_tensor_MO[i,A,j,B]
                exchange = eri_tensor_MO[i,j,A,B]
                
                H_singlet[ia,jb] = diag + 2*coulomb - exchange
                H_triplet[ia,jb] = diag - exchange

eig_s = LA.eigh(H_singlet)[0]
eig_t = LA.eigh(H_triplet)[0]

print("Singlet CIS:                    Triplet CIS:")
for n in range(min(10, len(eig_s))):
    print(f"  S{n+1}: {eig_s[n]:.6f} Ha ({eig_s[n]*27.2114:.2f} eV)"
          f"    T{n+1}: {eig_t[n]:.6f} Ha ({eig_t[n]*27.2114:.2f} eV)")

print("\n--- Starting Full CI Computation ---")
h_ao = mol.intor('int1e_kin') + mol.intor('int1e_nuc')  #'int1e_nuc' refers to interaction between N and e, not N and N
h_mo = C.T @ h_ao @ C
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]

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

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

def slater_codon(d1, d2):
    s1 = set(d1); s2 = set(d2)
    u1 = sorted(s1 - s2); u2 = sorted(s2 - s1)
    common = sorted(s1 & s2)
    def count_swaps(det, uniques):
        det = list(det)
        swaps = 0
        for k, u in enumerate(uniques):
            pos = det.index(u)
            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 = slater_codon(d1, d2)
    ndiff = len(u1)
    if ndiff == 0:
        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 H matrix...")
H_FCI = np.zeros((n_dets, n_dets))

for I in range(n_dets):
    for J in range(n_dets):
        val = H_element(dets[I], dets[J])
        H_FCI[I,J] = val
        H_FCI[J,I] = val

eigvals, eigvecs = LA.eigh(H_FCI)

print(f"FCI correlation: {eigvals[0] - E_elec:.8f} Hartree")

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_elec + enuc):.8f} Ha")