🔬 Chapter 3: Molecular Orbital Theory and Electronic Structure Calculations

import numpy as np import matplotlib.pyplot as plt from scipy.linalg import eigh class MinimalBasisH2: """SCF calculation for the H₂ molecule in a minimal basis (STO-1G)""" def __init__(self, R=1.4): """ R: internuclear distance (Bohr) """ self.R = R self.alpha = 1.0 # STO-1G exponent (simplified) def overlap_matrix(self): """Overlap matrix S""" R = self.R alpha = self.alpha # S_11 = S_22 = 1 (normalization) # S_12 = S_21 = ⟨φ_A|φ_B⟩ (simple approximation) S_12 = np.exp(-alpha * R) * (1 + alpha * R + (alpha * R)**2 / 3) S = np.array([[1.0, S_12], [S_12, 1.0]]) return S def core_hamiltonian(self): """Core Hamiltonian matrix H_core""" R = self.R alpha = self.alpha # Kinetic energy + nuclear attraction # H_11 = ⟨φ_A|T + V_A + V_B|φ_A⟩ H_11 = -alpha**2 / 2 - 1/R * (1 + 1/R) * np.exp(-2*alpha*R) - 1.0 # H_12 = ⟨φ_A|T + V_A + V_B|φ_B⟩ S_12 = np.exp(-alpha * R) * (1 + alpha * R + (alpha * R)**2 / 3) H_12 = -alpha**2 * S_12 / 2 - S_12 / R - S_12 / alpha H_core = np.array([[H_11, H_12], [H_12, H_11]]) return H_core def two_electron_integrals(self): """Two-electron integrals (simple approximation)""" R = self.R alpha = self.alpha # (μν|λσ) integrals (simplified) g = {} g['1111'] = 5/8 * alpha # ⟨φ_A φ_A|φ_A φ_A⟩ g['1122'] = 1/R * np.exp(-alpha * R) # ⟨φ_A φ_A|φ_B φ_B⟩ g['1212'] = 0.5 * g['1122'] # ⟨φ_A φ_B|φ_A φ_B⟩ g['2222'] = g['1111'] return g def build_fock_matrix(self, P, H_core, g): """Build the Fock matrix""" # F_μν = H_μν^core + Σ_λσ P_λσ [(μν|λσ) - 0.5(μλ|νσ)] F = H_core.copy() # Coulomb and exchange terms (simplified) F[0, 0] += P[0, 0] * g['1111'] + P[1, 1] * (g['1122'] - 0.5 * g['1212']) F[1, 1] += P[1, 1] * g['2222'] + P[0, 0] * (g['1122'] - 0.5 * g['1212']) F[0, 1] += P[0, 1] * (g['1212'] - 0.5 * g['1122']) F[1, 0] = F[0, 1] return F def scf_iteration(self, max_iter=50, conv_threshold=1e-6): """SCF iterative calculation""" S = self.overlap_matrix() H_core = self.core_hamiltonian() g = self.two_electron_integrals() # Initial density matrix (zero) P = np.zeros((2, 2)) energies = [] converged = False for iteration in range(max_iter): # Build the Fock matrix F = self.build_fock_matrix(P, H_core, g) # Solve the generalized eigenvalue problem: FC = SCε epsilon, C = eigh(F, S) # Update the density matrix (two-electron system, so one occupied orbital) P_new = 2 * np.outer(C[:, 0], C[:, 0]) # for two electrons # Electronic energy E_elec = 0.5 * np.sum(P_new * (H_core + F)) # Nuclear repulsion V_NN = 1 / self.R # Total energy E_total = E_elec + V_NN energies.append(E_total) # Convergence check if iteration > 0: delta_E = abs(E_total - energies[-2]) if delta_E < conv_threshold: converged = True break P = P_new return { 'converged': converged, 'iterations': iteration + 1, 'energy': E_total, 'orbital_energies': epsilon, 'coefficients': C, 'density_matrix': P, 'energy_history': energies } # SCF calculations at different internuclear distances R_range = np.linspace(0.5, 4.0, 30) energies_scf = [] orbital_energies_bonding = [] orbital_energies_antibonding = [] for R in R_range: h2 = MinimalBasisH2(R) result = h2.scf_iteration() energies_scf.append(result['energy']) orbital_energies_bonding.append(result['orbital_energies'][0]) orbital_energies_antibonding.append(result['orbital_energies'][1]) # Equilibrium internuclear distance E_array = np.array(energies_scf) R_eq_idx = np.argmin(E_array) R_eq = R_range[R_eq_idx] E_eq = E_array[R_eq_idx] # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Potential energy curve ax1 = axes[0, 0] ax1.plot(R_range, energies_scf, 'b-', linewidth=2, label='SCF Energy') ax1.plot(R_eq, E_eq, 'ro', markersize=10, label=f'R_eq = {R_eq:.2f} Bohr') ax1.axhline(0, color='k', linestyle='--', linewidth=1) ax1.set_xlabel('Internuclear distance R (Bohr)') ax1.set_ylabel('Energy (Hartree)') ax1.set_title('H₂ Potential Energy Curve (SCF)') ax1.legend() ax1.grid(True, alpha=0.3) # Orbital energies ax2 = axes[0, 1] ax2.plot(R_range, orbital_energies_bonding, 'b-', linewidth=2, label='Bonding MO') ax2.plot(R_range, orbital_energies_antibonding, 'r-', linewidth=2, label='Antibonding MO') ax2.axhline(0, color='k', linestyle='--', linewidth=1) ax2.set_xlabel('Internuclear distance R (Bohr)') ax2.set_ylabel('Orbital energy (Hartree)') ax2.set_title('Molecular Orbital Energies') ax2.legend() ax2.grid(True, alpha=0.3) # SCF convergence history (R = R_eq) ax3 = axes[1, 0] h2_eq = MinimalBasisH2(R_eq) result_eq = h2_eq.scf_iteration() ax3.plot(result_eq['energy_history'], 'go-', linewidth=2, markersize=6) ax3.set_xlabel('SCF iteration') ax3.set_ylabel('Total energy (Hartree)') ax3.set_title(f'SCF Convergence (R = {R_eq:.2f} Bohr)') ax3.grid(True, alpha=0.3) # Molecular orbital coefficients (R = R_eq) ax4 = axes[1, 1] C = result_eq['coefficients'] x = ['φ_A (H_A)', 'φ_B (H_B)'] width = 0.35 bonding_coeffs = C[:, 0] antibonding_coeffs = C[:, 1] x_pos = np.arange(len(x)) ax4.bar(x_pos - width/2, bonding_coeffs, width, label='Bonding', color='blue') ax4.bar(x_pos + width/2, antibonding_coeffs, width, label='Antibonding', color='red') ax4.set_ylabel('Coefficient') ax4.set_title(f'Molecular Orbital Coefficients (R = {R_eq:.2f} Bohr)') ax4.set_xticks(x_pos) ax4.set_xticklabels(x) ax4.legend() ax4.grid(True, alpha=0.3, axis='y') ax4.axhline(0, color='k', linewidth=0.5) plt.tight_layout() plt.savefig('qchem_scf_h2.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("=== SCF Calculation for the H₂ Molecule (Minimal Basis) ===\n") print(f"Equilibrium internuclear distance: R_eq = {R_eq:.3f} Bohr = {R_eq * 0.529:.3f} Å") print(f"Total energy: E = {E_eq:.6f} Hartree = {E_eq * 27.2114:.3f} eV") print(f"\nSCF convergence:") print(f" Iterations: {result_eq['iterations']}") print(f" Converged: {result_eq['converged']}") print(f"\nOrbital energies (R = R_eq):") print(f" Bonding MO: ε₁ = {result_eq['orbital_energies'][0]:.6f} Hartree") print(f" Antibonding MO: ε₂ = {result_eq['orbital_energies'][1]:.6f} Hartree")

import numpy as np import matplotlib.pyplot as plt def gaussian_1s(r, alpha): """Gaussian-type 1s orbital (normalized)""" N = (2 * alpha / np.pi)**(3/4) return N * np.exp(-alpha * r**2) def slater_1s(r, zeta=1.0): """Slater-type 1s orbital (normalized)""" N = (zeta**3 / np.pi)**0.5 return N * np.exp(-zeta * r) def sto_3g_1s(r): """STO-3G basis (approximating an STO with three Gaussians)""" # STO-3G parameters for the 1s orbital of the H atom (ζ=1.0) alphas = np.array([0.168856, 0.623913, 3.42525]) coeffs = np.array([0.444635, 0.535328, 0.154329]) result = np.zeros_like(r) for alpha, coeff in zip(alphas, coeffs): result += coeff * gaussian_1s(r, alpha) return result # Radial coordinate r = np.linspace(0, 5, 500) # Different basis functions psi_slater = slater_1s(r, zeta=1.0) psi_sto3g = sto_3g_1s(r) psi_gauss_single = gaussian_1s(r, alpha=0.3) # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Comparison of wavefunctions ax1 = axes[0, 0] ax1.plot(r, psi_slater, 'b-', linewidth=2, label='Slater 1s (ζ=1.0)') ax1.plot(r, psi_sto3g, 'r--', linewidth=2, label='STO-3G') ax1.plot(r, psi_gauss_single, 'g:', linewidth=2, label='Single Gaussian (α=0.3)') ax1.set_xlabel('r (Bohr)') ax1.set_ylabel('ψ(r)') ax1.set_title('Comparison of Basis Functions') ax1.legend() ax1.grid(True, alpha=0.3) # Radial probability density ax2 = axes[0, 1] ax2.plot(r, r**2 * psi_slater**2, 'b-', linewidth=2, label='Slater') ax2.plot(r, r**2 * psi_sto3g**2, 'r--', linewidth=2, label='STO-3G') ax2.set_xlabel('r (Bohr)') ax2.set_ylabel('r² |ψ(r)|²') ax2.set_title('Radial Probability Density') ax2.legend() ax2.grid(True, alpha=0.3) # Constituent Gaussians of STO-3G ax3 = axes[1, 0] alphas_sto3g = np.array([0.168856, 0.623913, 3.42525]) coeffs_sto3g = np.array([0.444635, 0.535328, 0.154329]) for i, (alpha, coeff) in enumerate(zip(alphas_sto3g, coeffs_sto3g)): psi_component = coeff * gaussian_1s(r, alpha) ax3.plot(r, psi_component, linewidth=2, label=f'G{i+1} (α={alpha:.2f}, c={coeff:.3f})') ax3.plot(r, psi_sto3g, 'k-', linewidth=3, label='STO-3G (sum)') ax3.set_xlabel('r (Bohr)') ax3.set_ylabel('ψ(r)') ax3.set_title('Constituent Gaussians of STO-3G') ax3.legend() ax3.grid(True, alpha=0.3) # Effect of basis set size (conceptual) ax4 = axes[1, 1] basis_sets = ['Minimal\n(STO-3G)', 'Double-ζ\n(6-31G)', 'Triple-ζ\n(6-311G)', 'cc-pVDZ', 'cc-pVTZ'] n_functions = [1, 2, 3, 5, 14] # Number of functions for the H atom (approximate) relative_accuracy = [0.8, 0.92, 0.96, 0.98, 0.995] x_pos = np.arange(len(basis_sets)) bars = ax4.bar(x_pos, relative_accuracy, color=['red', 'orange', 'yellow', 'lightgreen', 'green']) ax4.set_ylabel('Relative accuracy') ax4.set_title('Basis Set Size and Accuracy') ax4.set_xticks(x_pos) ax4.set_xticklabels(basis_sets, rotation=15, ha='right') ax4.set_ylim([0.7, 1.0]) ax4.grid(True, alpha=0.3, axis='y') for i, (bar, n_func) in enumerate(zip(bars, n_functions)): height = bar.get_height() ax4.text(bar.get_x() + bar.get_width()/2., height + 0.01, f'{n_func} funcs', ha='center', va='bottom', fontsize=9) plt.tight_layout() plt.savefig('qchem_basis_functions.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("\n=== Basis Functions ===\n") print("STO-3G parameters (H 1s):") for i, (alpha, coeff) in enumerate(zip(alphas_sto3g, coeffs_sto3g)): print(f" Gaussian {i+1}: α = {alpha:.6f}, c = {coeff:.6f}") print("\nGuidelines for choosing a basis set:") print(" - Minimal (STO-3G): qualitative understanding, large systems") print(" - Double-ζ (6-31G): standard calculations, reasonable accuracy") print(" - Triple-ζ (6-311G): high-accuracy calculations") print(" - cc-pVnZ: benchmark calculations, correlation energy")

import numpy as np import matplotlib.pyplot as plt from scipy.linalg import eigh class HuckelMO: """Hückel molecular orbital method""" def __init__(self, adjacency_matrix, alpha=0, beta=-1): """ adjacency_matrix: adjacency matrix (presence/absence of bonds) alpha: Coulomb integral (energy zero point) beta: resonance integral """ self.adjacency = adjacency_matrix self.n_atoms = len(adjacency_matrix) self.alpha = alpha self.beta = beta def build_hamiltonian(self): """Hückel Hamiltonian matrix""" H = self.alpha * np.eye(self.n_atoms) + self.beta * self.adjacency return H def solve(self): """Solve the eigenvalue problem""" H = self.build_hamiltonian() eigenvalues, eigenvectors = eigh(H) # Sort by energy idx = np.argsort(eigenvalues) eigenvalues = eigenvalues[idx] eigenvectors = eigenvectors[:, idx] return eigenvalues, eigenvectors def pi_energy(self, n_electrons): """π-electron energy""" eigenvalues, _ = self.solve() # Sum of the energies of the occupied orbitals n_occupied = n_electrons // 2 # two electrons per orbital E_pi = 2 * np.sum(eigenvalues[:n_occupied]) return E_pi # Adjacency matrix of benzene (C6H6) benzene_adj = np.array([ [0, 1, 0, 0, 0, 1], [1, 0, 1, 0, 0, 0], [0, 1, 0, 1, 0, 0], [0, 0, 1, 0, 1, 0], [0, 0, 0, 1, 0, 1], [1, 0, 0, 0, 1, 0] ]) # Adjacency matrix of butadiene (C4H6) butadiene_adj = np.array([ [0, 1, 0, 0], [1, 0, 1, 0], [0, 1, 0, 1], [0, 0, 1, 0] ]) # Ethylene (C2H4) ethylene_adj = np.array([ [0, 1], [1, 0] ]) # Calculation for benzene benzene = HuckelMO(benzene_adj, alpha=0, beta=-1) E_benzene, C_benzene = benzene.solve() E_pi_benzene = benzene.pi_energy(6) # 6 π electrons # Calculation for butadiene butadiene = HuckelMO(butadiene_adj, alpha=0, beta=-1) E_butadiene, C_butadiene = butadiene.solve() # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Energy level diagram for benzene ax1 = axes[0, 0] colors_benzene = ['blue', 'green', 'green', 'red', 'red', 'red'] occupancy_benzene = [2, 2, 2, 0, 0, 0] # number of electrons in each level for i, (E, occ, color) in enumerate(zip(E_benzene, occupancy_benzene, colors_benzene)): ax1.hlines(E, i-0.3, i+0.3, colors=color, linewidth=3) # Show the electrons if occ > 0: ax1.plot([i-0.1, i+0.1], [E, E], 'o', color='black', markersize=8) ax1.axhline(0, color='k', linestyle='--', linewidth=1, label='α (reference)') ax1.set_xticks(range(len(E_benzene))) ax1.set_xticklabels([f'MO{i+1}' for i in range(len(E_benzene))]) ax1.set_ylabel('Energy (units of β)') ax1.set_title('Hückel MO Levels of Benzene') ax1.legend() ax1.grid(True, alpha=0.3, axis='y') # Energy levels of butadiene ax2 = axes[0, 1] for i, E in enumerate(E_butadiene): color = 'blue' if i < 2 else 'red' ax2.hlines(E, i-0.3, i+0.3, colors=color, linewidth=3) if i < 2: # occupied orbitals ax2.plot([i-0.1, i+0.1], [E, E], 'o', color='black', markersize=8) ax2.axhline(0, color='k', linestyle='--', linewidth=1) ax2.set_xticks(range(len(E_butadiene))) ax2.set_xticklabels([f'π{i+1}' for i in range(len(E_butadiene))]) ax2.set_ylabel('Energy (units of β)') ax2.set_title('Hückel MO Levels of Butadiene') ax2.grid(True, alpha=0.3, axis='y') # Molecular orbital coefficients of benzene (HOMO) ax3 = axes[1, 0] homo_index = 2 # 3rd MO (0-indexed) homo_coeffs = C_benzene[:, homo_index] theta = np.linspace(0, 2*np.pi, 7) x_pos = np.cos(theta[:6]) y_pos = np.sin(theta[:6]) # Color-code by the sign of the coefficient colors_coeff = ['red' if c > 0 else 'blue' for c in homo_coeffs] sizes = np.abs(homo_coeffs) * 500 ax3.scatter(x_pos, y_pos, s=sizes, c=colors_coeff, alpha=0.6, edgecolors='black', linewidth=2) # Draw the benzene ring for i in range(6): ax3.plot([x_pos[i], x_pos[(i+1)%6]], [y_pos[i], y_pos[(i+1)%6]], 'k-', linewidth=1) ax3.set_xlim([-1.5, 1.5]) ax3.set_ylim([-1.5, 1.5]) ax3.set_aspect('equal') ax3.set_title(f'Benzene HOMO (MO{homo_index+1}) Coefficients') ax3.set_xticks([]) ax3.set_yticks([]) ax3.text(0, -1.8, 'Red: positive, Blue: negative', ha='center', fontsize=10) # Dependence of π-electron energy on conjugated chain length ax4 = axes[1, 1] chain_lengths = range(2, 11) pi_energies = [] for n in chain_lengths: # Adjacency matrix of a linear conjugated system adj = np.diag(np.ones(n-1), 1) + np.diag(np.ones(n-1), -1) mol = HuckelMO(adj, alpha=0, beta=-1) E_pi = mol.pi_energy(n) # n π electrons pi_energies.append(E_pi) ax4.plot(chain_lengths, pi_energies, 'go-', linewidth=2, markersize=8) ax4.set_xlabel('Number of carbon atoms') ax4.set_ylabel('Total π energy (units of β)') ax4.set_title('Conjugated Chain Length and π-Electron Energy') ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('qchem_huckel_mo.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("\n=== Hückel Molecular Orbital Method ===\n") print("Benzene (C6H6):") print(f" MO energies (α + xβ):") for i, E in enumerate(E_benzene): print(f" MO{i+1}: E = α + ({E:.4f})β") print(f" Total π-electron energy: E_π = {E_pi_benzene:.4f}β") print(f" Resonance stabilization energy: {E_pi_benzene - 6*(-1):.4f}β") print("\nButadiene (C4H6):") for i, E in enumerate(E_butadiene): print(f" π{i+1}: E = α + ({E:.4f})β")

import numpy as np import matplotlib.pyplot as plt from scipy.linalg import eigh class SimpleCIS: """Simplified CIS (Configuration Interaction Singles)""" def __init__(self, n_orbitals=4): """ n_orbitals: number of molecular orbitals """ self.n_orbitals = n_orbitals # Dummy orbital energies (in Hartree) self.orbital_energies = np.array([ -0.5, # HOMO-1 -0.3, # HOMO 0.2, # LUMO 0.4 # LUMO+1 ]) # Number of electrons (assume a two-electron system) self.n_electrons = 2 self.homo_index = 1 # 0-indexed def ground_state_energy(self): """Ground-state energy (HF)""" # Sum of the energies of the occupied orbitals E_HF = 2 * np.sum(self.orbital_energies[:self.n_electrons//2]) return E_HF def single_excitations(self): """Single-electron excited configurations""" excitations = [] # HOMO → LUMO, HOMO → LUMO+1, etc. for occ in range(self.n_electrons // 2): for virt in range(self.n_electrons // 2, self.n_orbitals): excitation_energy = self.orbital_energies[virt] - self.orbital_energies[occ] excitations.append({ 'from': occ, 'to': virt, 'energy': excitation_energy }) return excitations def cis_matrix(self): """CIS Hamiltonian matrix (simplified version)""" excitations = self.single_excitations() n_exc = len(excitations) H_CIS = np.zeros((n_exc, n_exc)) # Diagonal elements: excitation energies for i, exc in enumerate(excitations): H_CIS[i, i] = exc['energy'] # Off-diagonal elements: configuration interaction (simplified: neglected) # In practice, Coulomb and exchange integrals are computed return H_CIS, excitations def solve_cis(self): """Solve the CIS equations""" H_CIS, excitations = self.cis_matrix() eigenvalues, eigenvectors = eigh(H_CIS) return eigenvalues, eigenvectors, excitations # CIS calculation cis = SimpleCIS(n_orbitals=4) E_ground = cis.ground_state_energy() exc_energies, exc_states, excitations = cis.solve_cis() # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Orbital energy diagram ax1 = axes[0, 0] orb_energies = cis.orbital_energies colors_orb = ['blue', 'blue', 'red', 'red'] labels_orb = ['HOMO-1', 'HOMO', 'LUMO', 'LUMO+1'] for i, (E, color, label) in enumerate(zip(orb_energies, colors_orb, labels_orb)): ax1.hlines(E, i-0.3, i+0.3, colors=color, linewidth=3) if i < 2: # occupied orbitals ax1.plot([i-0.1, i+0.1], [E, E], 'o', color='black', markersize=8) ax1.text(i, E - 0.15, label, ha='center', fontsize=9) ax1.axhline(0, color='k', linestyle='--', linewidth=1, label='Vacuum level') ax1.set_xlim([-0.5, 3.5]) ax1.set_ylabel('Energy (Hartree)') ax1.set_title('Molecular Orbital Energy Levels') ax1.set_xticks([]) ax1.legend() ax1.grid(True, alpha=0.3, axis='y') # Excited-state energies ax2 = axes[0, 1] state_labels = [f"S{i+1}" for i in range(len(exc_energies))] ax2.barh(range(len(exc_energies)), exc_energies, color='orange', edgecolor='black') for i, (E, label) in enumerate(zip(exc_energies, state_labels)): ax2.text(E + 0.02, i, f'{label}: {E:.3f} Ha', va='center', fontsize=10) ax2.set_xlabel('Excitation energy (Hartree)') ax2.set_ylabel('Excited state') ax2.set_title('CIS Excited States') ax2.set_yticks(range(len(exc_energies))) ax2.set_yticklabels(state_labels) ax2.grid(True, alpha=0.3, axis='x') # Composition of an excited configuration (first excited state) ax3 = axes[1, 0] state_idx = 0 state_vector = exc_states[:, state_idx] x_pos = np.arange(len(state_vector)) ax3.bar(x_pos, np.abs(state_vector)**2, color='green', edgecolor='black') exc_labels = [f"{exc['from']}→{exc['to']}" for exc in excitations] ax3.set_xticks(x_pos) ax3.set_xticklabels(exc_labels) ax3.set_xlabel('Excitation') ax3.set_ylabel('|Coefficient|²') ax3.set_title(f'Composition of Excited State S{state_idx+1}') ax3.grid(True, alpha=0.3, axis='y') # Conceptual diagram of correlation energy ax4 = axes[1, 1] methods = ['HF', 'CIS', 'CISD', 'CCSD', 'FCI'] relative_corr = [0, 0.1, 0.5, 0.8, 1.0] # relative correlation energy recovery colors_method = ['red', 'orange', 'yellow', 'lightgreen', 'green'] bars = ax4.barh(range(len(methods)), relative_corr, color=colors_method, edgecolor='black') ax4.set_xlabel('Correlation energy recovery') ax4.set_ylabel('Method') ax4.set_title('Treatment of Electron Correlation (Conceptual)') ax4.set_yticks(range(len(methods))) ax4.set_yticklabels(methods) ax4.set_xlim([0, 1.1]) ax4.grid(True, alpha=0.3, axis='x') for i, (bar, corr) in enumerate(zip(bars, relative_corr)): width = bar.get_width() ax4.text(width + 0.02, bar.get_y() + bar.get_height()/2, f'{corr*100:.0f}%', va='center', fontsize=10) plt.tight_layout() plt.savefig('qchem_electron_correlation.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("\n=== Electron Correlation and CI ===\n") print(f"Ground-state energy (HF): {E_ground:.6f} Hartree") print(f"\nExcited states (CIS):") for i, E_exc in enumerate(exc_energies): print(f" S{i+1}: ΔE = {E_exc:.6f} Hartree = {E_exc * 27.2114:.3f} eV") print("\nHierarchy of post-HF methods:") print(" HF: mean-field approximation, no electron correlation") print(" CIS: single excitations only, excited-state calculation") print(" CISD: single and double excitations, part of the correlation") print(" CCSD: Coupled Cluster, high accuracy") print(" FCI: Full CI, exact solution (small molecules only)")