Chapter 5: Practical Quantum Chemistry Calculations

import numpy as np import matplotlib.pyplot as plt from scipy.linalg import eigh # Minimal PySCF-compatible implementation (for demonstration) class WaterMoleculeHF: """Hartree-Fock calculation of the water molecule (simplified version)""" def __init__(self, basis='sto-3g'): # Water molecule geometry (O-H distance 0.96 Å, H-O-H angle 104.5°) self.atoms = [ ('O', [0.0, 0.0, 0.0]), ('H', [0.0, 0.757, 0.587]), ('H', [0.0, -0.757, 0.587]) ] self.n_electrons = 10 # 8(O) + 1(H) + 1(H) self.n_basis = 7 # STO-3G: O(5) + H(1) + H(1) def compute_integrals(self): """Compute the integral matrices (simplified version)""" # Overlap integral matrix (nearly identity after diagonalization) S = np.eye(self.n_basis) for i in range(self.n_basis): for j in range(i+1, self.n_basis): S[i,j] = S[j,i] = 0.2 * np.exp(-0.5 * (i-j)**2) # Core Hamiltonian (kinetic energy + nuclear-electron attraction) H_core = -np.diag([20.0, 15.0, 15.0, 15.0, 15.0, 7.0, 7.0]) # Two-electron integrals (simplified version) g = np.zeros((self.n_basis, self.n_basis, self.n_basis, self.n_basis)) for i in range(self.n_basis): for j in range(self.n_basis): g[i,i,j,j] = 10.0 / (i+j+2) g[i,j,i,j] = 8.0 / (i+j+3) return S, H_core, g def build_fock(self, P, H_core, g): """Construct the Fock matrix""" F = H_core.copy() for i in range(self.n_basis): for j in range(self.n_basis): for k in range(self.n_basis): for l in range(self.n_basis): F[i,j] += P[k,l] * (g[i,j,k,l] - 0.5 * g[i,k,j,l]) return F def run_scf(self, max_iter=30, conv=1e-6): """SCF calculation""" S, H_core, g = self.compute_integrals() # Initial density matrix P = np.zeros((self.n_basis, self.n_basis)) E_old = 0.0 energies = [] for iteration in range(max_iter): # Build the Fock matrix F = self.build_fock(P, H_core, g) # Generalized eigenvalue problem epsilon, C = eigh(F, S) # Update the density matrix (10 electrons = 5 occupied orbitals) n_occ = self.n_electrons // 2 P_new = 2 * C[:, :n_occ] @ C[:, :n_occ].T # Energy calculation E_elec = np.sum(P_new * (H_core + F)) / 2 E_nuc = 9.0 # Nuclear repulsion (simplified value) E_total = E_elec + E_nuc energies.append(E_total) # Convergence check if np.abs(E_total - E_old) < conv: print(f"SCF converged: {iteration+1} iterations") break E_old = E_total P = P_new return E_total, epsilon, C, energies # Run the HF calculation for the water molecule h2o = WaterMoleculeHF() E_total, epsilon, C, energies = h2o.run_scf() # Visualization fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # (1) SCF convergence axes[0,0].plot(energies, 'o-', linewidth=2) axes[0,0].set_xlabel('Iteration', fontsize=12) axes[0,0].set_ylabel('Energy (Hartree)', fontsize=12) axes[0,0].set_title('SCF Convergence', fontsize=14, fontweight='bold') axes[0,0].grid(True, alpha=0.3) # (2) Orbital energy levels n_occ = h2o.n_electrons // 2 colors = ['blue' if i < n_occ else 'red' for i in range(h2o.n_basis)] axes[0,1].barh(range(h2o.n_basis), epsilon, color=colors, alpha=0.7) axes[0,1].axvline(0, color='black', linestyle='--', alpha=0.5) axes[0,1].set_xlabel('Energy (Hartree)', fontsize=12) axes[0,1].set_ylabel('Orbital Index', fontsize=12) axes[0,1].set_title('Orbital Energy Levels', fontsize=14, fontweight='bold') axes[0,1].legend(['Occupied', 'Virtual'], loc='best') # (3) HOMO-LUMO coefficients homo_idx = n_occ - 1 lumo_idx = n_occ x = np.arange(h2o.n_basis) axes[1,0].bar(x - 0.2, C[:, homo_idx], width=0.4, label='HOMO', alpha=0.8) axes[1,0].bar(x + 0.2, C[:, lumo_idx], width=0.4, label='LUMO', alpha=0.8) axes[1,0].set_xlabel('Basis Function', fontsize=12) axes[1,0].set_ylabel('Coefficient', fontsize=12) axes[1,0].set_title('HOMO-LUMO Molecular Orbitals', fontsize=14, fontweight='bold') axes[1,0].legend() axes[1,0].grid(True, alpha=0.3, axis='y') # (4) Summary of results axes[1,1].axis('off') summary = f""" H₂O Hartree-Fock Calculation Total Energy: {E_total:.6f} Hartree {E_total * 27.211:.2f} eV HOMO Energy: {epsilon[homo_idx]:.6f} Hartree LUMO Energy: {epsilon[lumo_idx]:.6f} Hartree HOMO-LUMO Gap: {(epsilon[lumo_idx] - epsilon[homo_idx]) * 27.211:.2f} eV Basis Set: STO-3G Electrons: {h2o.n_electrons} Basis Funcs: {h2o.n_basis} """ axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace', verticalalignment='center') plt.tight_layout() plt.savefig('water_hf_calculation.png', dpi=150, bbox_inches='tight') plt.show() print(f"H₂O Total Energy: {E_total:.6f} Hartree") print(f"HOMO-LUMO Gap: {(epsilon[lumo_idx] - epsilon[homo_idx]) * 27.211:.2f} eV")

class WaterMoleculeDFT(WaterMoleculeHF): """DFT calculation of the water molecule (simplified B3LYP version)""" def __init__(self, functional='b3lyp'): super().__init__() self.functional = functional def xc_energy(self, rho): """Exchange-correlation energy (B3LYP approximation)""" # LDA exchange C_x = (3/4) * (3/np.pi)**(1/3) E_x_lda = -C_x * rho**(4/3) # GGA correction (gradient approximation) E_x_gga = E_x_lda * (1.0 + 0.1 * rho**(1/3)) # B3LYP mixing (20% HF + 80% DFT) E_xc = 0.8 * E_x_gga return E_xc def build_fock_dft(self, P, H_core, g): """Kohn-Sham Fock matrix""" F = H_core.copy() # Coulomb + Exchange for i in range(self.n_basis): for j in range(self.n_basis): for k in range(self.n_basis): for l in range(self.n_basis): # Coulomb F[i,j] += P[k,l] * g[i,j,k,l] # Exchange (20% HF for B3LYP) F[i,j] -= 0.2 * P[k,l] * g[i,k,j,l] # XC potential (simplified) rho = np.diag(P) V_xc = self.xc_energy(rho + 1e-10) F += np.diag(V_xc) return F def run_dft(self, max_iter=30, conv=1e-6): """DFT-SCF calculation""" S, H_core, g = self.compute_integrals() P = np.zeros((self.n_basis, self.n_basis)) E_old = 0.0 energies = [] for iteration in range(max_iter): F = self.build_fock_dft(P, H_core, g) epsilon, C = eigh(F, S) n_occ = self.n_electrons // 2 P_new = 2 * C[:, :n_occ] @ C[:, :n_occ].T E_elec = np.sum(P_new * (H_core + F)) / 2 E_nuc = 9.0 E_total = E_elec + E_nuc energies.append(E_total) if np.abs(E_total - E_old) < conv: print(f"DFT-SCF converged: {iteration+1} iterations") break E_old = E_total P = P_new return E_total, epsilon, C, energies # Run the DFT calculation h2o_dft = WaterMoleculeDFT() E_dft, eps_dft, C_dft, energies_dft = h2o_dft.run_dft() # HF vs DFT comparison print(f"\n=== H₂O Energy Comparison ===") print(f"HF Energy: {E_total:.6f} Hartree") print(f"DFT Energy: {E_dft:.6f} Hartree") print(f"Difference: {(E_dft - E_total)*1000:.2f} mHartree") n_occ = h2o.n_electrons // 2 homo_hf = epsilon[n_occ-1] lumo_hf = epsilon[n_occ] homo_dft = eps_dft[n_occ-1] lumo_dft = eps_dft[n_occ] print(f"\nHOMO-LUMO Gap:") print(f"HF: {(lumo_hf - homo_hf)*27.211:.2f} eV") print(f"DFT: {(lumo_dft - homo_dft)*27.211:.2f} eV")

from scipy.optimize import minimize def h2_energy(R): """Energy of the H₂ molecule (Morse potential approximation)""" D_e = 4.75 # eV alpha = 1.44 # Å⁻¹ R_e = 0.741 # Å E = D_e * ((1 - np.exp(-alpha * (R - R_e)))**2 - 1) return E def h2_gradient(R): """Energy gradient""" D_e = 4.75 alpha = 1.44 R_e = 0.741 exp_term = np.exp(-alpha * (R - R_e)) dE_dR = 2 * D_e * alpha * (1 - exp_term) * exp_term return dE_dR # Geometry optimization (initial value: R=1.0 Å) result = minimize(h2_energy, x0=[1.0], method='BFGS', jac=h2_gradient) R_opt = result.x[0] E_opt = result.fun print(f"Optimization results:") print(f" Equilibrium distance: {R_opt:.4f} Å") print(f" Bond energy: {E_opt:.4f} eV") # Potential curve and PES R_range = np.linspace(0.3, 3.0, 100) E_range = [h2_energy(r) for r in R_range] fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # (1) Potential energy curve axes[0,0].plot(R_range, E_range, 'b-', linewidth=2, label='Morse Potential') axes[0,0].plot(R_opt, E_opt, 'ro', markersize=10, label=f'Minimum: {R_opt:.3f}Å') axes[0,0].axhline(0, color='k', linestyle='--', alpha=0.3) axes[0,0].set_xlabel('H-H Distance (Å)', fontsize=12) axes[0,0].set_ylabel('Energy (eV)', fontsize=12) axes[0,0].set_title('H₂ Potential Energy Curve', fontsize=14, fontweight='bold') axes[0,0].legend() axes[0,0].grid(True, alpha=0.3) # (2) Energy gradient gradient_range = [h2_gradient(r) for r in R_range] axes[0,1].plot(R_range, gradient_range, 'g-', linewidth=2) axes[0,1].axhline(0, color='r', linestyle='--', linewidth=2, label='Zero gradient') axes[0,1].plot(R_opt, 0, 'ro', markersize=10) axes[0,1].set_xlabel('H-H Distance (Å)', fontsize=12) axes[0,1].set_ylabel('dE/dR (eV/Å)', fontsize=12) axes[0,1].set_title('Energy Gradient', fontsize=14, fontweight='bold') axes[0,1].legend() axes[0,1].grid(True, alpha=0.3) # (3) Vibrational energy levels # Morse oscillator levels omega = 4401.2 # cm⁻¹ (H₂ vibrational frequency) D_e_cm = D_e * 8065.5 # eV → cm⁻¹ n_max = int(D_e_cm / omega) n_levels = np.arange(0, min(n_max, 15)) E_vib = omega * (n_levels + 0.5) - omega**2 * (n_levels + 0.5)**2 / (4 * D_e_cm) axes[1,0].barh(n_levels, E_vib / 8065.5, height=0.8, alpha=0.7) axes[1,0].set_ylabel('Vibrational Quantum Number (n)', fontsize=12) axes[1,0].set_xlabel('Energy (eV)', fontsize=12) axes[1,0].set_title('Vibrational Energy Levels', fontsize=14, fontweight='bold') axes[1,0].grid(True, alpha=0.3, axis='x') # (4) Optimization summary axes[1,1].axis('off') summary = f""" H₂ Geometry Optimization Equilibrium Distance: {R_opt:.4f} Å Bond Energy: {-E_opt:.4f} eV {-E_opt * 96.485:.1f} kJ/mol Vibrational Frequency: {omega:.1f} cm⁻¹ Zero-Point Energy: {E_vib[0] / 8065.5:.4f} eV Optimization Method: BFGS Convergence: {result.success} Iterations: {result.nit} """ axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace', verticalalignment='center') plt.tight_layout() plt.savefig('h2_geometry_optimization.png', dpi=150, bbox_inches='tight') plt.show()

class EthyleneTDDFT: """TD-DFT calculation of the ethylene molecule (simplified version)""" def __init__(self): self.n_occ = 6 # Number of occupied orbitals self.n_virt = 4 # Number of virtual orbitals # Orbital energies (approximated from experiment, in Hartree units) self.epsilon_occ = np.array([-0.8, -0.7, -0.6, -0.5, -0.4, -0.35]) self.epsilon_virt = np.array([0.1, 0.2, 0.3, 0.4]) def compute_excitations(self): """One-electron excitation energies""" excitations = [] for i, e_occ in enumerate(self.epsilon_occ): for a, e_virt in enumerate(self.epsilon_virt): omega = e_virt - e_occ # Oscillator strength (simplified approximation) f = 0.1 * np.exp(-2 * (i + a)) excitations.append({ 'from': i, 'to': a, 'energy_hartree': omega, 'energy_eV': omega * 27.211, 'wavelength_nm': 1239.8 / (omega * 27.211), 'oscillator_strength': f }) return sorted(excitations, key=lambda x: x['energy_eV']) def absorption_spectrum(self, excitations, lambda_range): """Compute the absorption spectrum""" spectrum = np.zeros_like(lambda_range) for exc in excitations: # Lorentzian broadening gamma = 20 # nm (full width at half maximum) lambda_0 = exc['wavelength_nm'] f = exc['oscillator_strength'] spectrum += f * (gamma / 2) / ((lambda_range - lambda_0)**2 + (gamma / 2)**2) return spectrum # Run the TD-DFT calculation ethylene = EthyleneTDDFT() excitations = ethylene.compute_excitations() # Major excited states (oscillator strength > 0.01) major_exc = [exc for exc in excitations if exc['oscillator_strength'] > 0.01] print("Major excited states:") for i, exc in enumerate(major_exc[:5], 1): print(f"{i}. {exc['energy_eV']:.2f} eV ({exc['wavelength_nm']:.1f} nm) " f"[{exc['from']}→{exc['to']}] f={exc['oscillator_strength']:.3f}") # Compute the spectrum lambda_range = np.linspace(100, 400, 500) spectrum = ethylene.absorption_spectrum(excitations, lambda_range) # Visualization fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # (1) Orbital energy diagram y_occ = np.arange(ethylene.n_occ) y_virt = np.arange(ethylene.n_virt) + ethylene.n_occ + 1 axes[0,0].barh(y_occ, ethylene.epsilon_occ * 27.211, height=0.6, color='blue', alpha=0.7, label='Occupied') axes[0,0].barh(y_virt, ethylene.epsilon_virt * 27.211, height=0.6, color='red', alpha=0.7, label='Virtual') axes[0,0].axvline(0, color='k', linestyle='--', alpha=0.3) # Arrow for the π→π* transition homo = ethylene.n_occ - 1 lumo = ethylene.n_occ + 1 axes[0,0].annotate('', xy=(ethylene.epsilon_virt[0]*27.211, lumo), xytext=(ethylene.epsilon_occ[-1]*27.211, homo), arrowprops=dict(arrowstyle='->', lw=2, color='green')) axes[0,0].text(-5, (homo+lumo)/2, 'π→π*', fontsize=12, color='green', fontweight='bold') axes[0,0].set_xlabel('Energy (eV)', fontsize=12) axes[0,0].set_ylabel('Orbital Index', fontsize=12) axes[0,0].set_title('Molecular Orbital Diagram', fontsize=14, fontweight='bold') axes[0,0].legend() # (2) Excitation energy spectrum (stick spectrum) energies = [exc['energy_eV'] for exc in major_exc] intensities = [exc['oscillator_strength'] for exc in major_exc] axes[0,1].stem(energies, intensities, basefmt=' ') axes[0,1].set_xlabel('Excitation Energy (eV)', fontsize=12) axes[0,1].set_ylabel('Oscillator Strength', fontsize=12) axes[0,1].set_title('Excitation Spectrum', fontsize=14, fontweight='bold') axes[0,1].grid(True, alpha=0.3) # (3) UV-Vis absorption spectrum axes[1,0].plot(lambda_range, spectrum, 'b-', linewidth=2) axes[1,0].fill_between(lambda_range, spectrum, alpha=0.3) axes[1,0].set_xlabel('Wavelength (nm)', fontsize=12) axes[1,0].set_ylabel('Absorbance (a.u.)', fontsize=12) axes[1,0].set_title('UV-Vis Absorption Spectrum', fontsize=14, fontweight='bold') axes[1,0].set_xlim(100, 400) axes[1,0].grid(True, alpha=0.3) # Annotate the position of the main peak max_idx = np.argmax(spectrum) max_lambda = lambda_range[max_idx] axes[1,0].annotate(f'{max_lambda:.1f} nm', xy=(max_lambda, spectrum[max_idx]), xytext=(max_lambda+50, spectrum[max_idx]*0.8), arrowprops=dict(arrowstyle='->', color='red'), fontsize=11, color='red', fontweight='bold') # (4) Excited-state summary axes[1,1].axis('off') summary = f""" C₂H₄ TD-DFT Excitation Analysis Lowest Excitation: Energy: {major_exc[0]['energy_eV']:.2f} eV Wavelength: {major_exc[0]['wavelength_nm']:.1f} nm Type: π → π* f: {major_exc[0]['oscillator_strength']:.3f} Max Absorption: λmax: {max_lambda:.1f} nm Region: {'UV-C' if max_lambda < 280 else 'UV-B' if max_lambda < 315 else 'UV-A'} Total Excitations: {len(excitations)} Major (f>0.01): {len(major_exc)} """ axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace', verticalalignment='center') plt.tight_layout() plt.savefig('ethylene_tddft.png', dpi=150, bbox_inches='tight') plt.show()

class KronigPenney1D: """One-dimensional periodic potential (Kronig-Penney model)""" def __init__(self, a=1.0, V0=10.0, N_pw=21): self.a = a # Lattice constant self.V0 = V0 # Potential depth self.N_pw = N_pw # Number of plane waves def potential(self, x): """Periodic potential (cosine type)""" return -self.V0 * np.cos(2 * np.pi * x / self.a) def hamiltonian_matrix(self, k): """Hamiltonian matrix at a k point""" # Plane-wave basis: exp(i(k + n*G)x), G = 2π/a G = 2 * np.pi / self.a n_indices = np.arange(-self.N_pw // 2, self.N_pw // 2 + 1) H = np.zeros((self.N_pw, self.N_pw), dtype=complex) for i, n_i in enumerate(n_indices): for j, n_j in enumerate(n_indices): k_i = k + n_i * G k_j = k + n_j * G if i == j: # Kinetic energy H[i,j] = 0.5 * k_i**2 elif abs(n_i - n_j) == 1: # Potential (first-nearest neighbor only) H[i,j] = -self.V0 / 2 return H def compute_bands(self, k_points): """Band structure calculation""" bands = [] for k in k_points: H = self.hamiltonian_matrix(k) eigenvalues = eigh(H, eigvals_only=True) bands.append(eigenvalues) return np.array(bands) # Band structure calculation kp = KronigPenney1D(a=1.0, V0=5.0, N_pw=21) # First Brillouin zone: -π/a ≤ k ≤ π/a k_points = np.linspace(-np.pi/kp.a, np.pi/kp.a, 100) bands = kp.compute_bands(k_points) # Band gap detection def find_band_gap(bands): """Detect the direct band gap""" valence_band_max = np.max(bands[:, :11]) # Lower bands conduction_band_min = np.min(bands[:, 11:]) # Upper bands gap = conduction_band_min - valence_band_max return gap, valence_band_max, conduction_band_min gap, vbm, cbm = find_band_gap(bands) # Visualization fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # (1) Periodic potential x = np.linspace(0, 3*kp.a, 500) V = kp.potential(x) axes[0,0].plot(x, V, 'b-', linewidth=2) for i in range(4): axes[0,0].axvline(i*kp.a, color='gray', linestyle='--', alpha=0.5) axes[0,0].set_xlabel('Position x/a', fontsize=12) axes[0,0].set_ylabel('Potential V(x)', fontsize=12) axes[0,0].set_title('Periodic Potential', fontsize=14, fontweight='bold') axes[0,0].grid(True, alpha=0.3) # (2) Band structure for n in range(min(10, bands.shape[1])): color = 'blue' if n < 5 else 'red' axes[0,1].plot(k_points * kp.a / np.pi, bands[:, n], color=color, linewidth=2) axes[0,1].axhline(vbm, color='blue', linestyle='--', alpha=0.5, label='VBM') axes[0,1].axhline(cbm, color='red', linestyle='--', alpha=0.5, label='CBM') axes[0,1].axvline(0, color='black', linestyle='-', alpha=0.3) axes[0,1].set_xlabel('Wave Vector (k/π/a)', fontsize=12) axes[0,1].set_ylabel('Energy (a.u.)', fontsize=12) axes[0,1].set_title('Band Structure', fontsize=14, fontweight='bold') axes[0,1].set_xlim(-1, 1) axes[0,1].legend() axes[0,1].grid(True, alpha=0.3) # Highlight the band gap region axes[0,1].fill_between([-1, 1], vbm, cbm, alpha=0.2, color='yellow') # (3) Density of states (DOS) energy_bins = np.linspace(bands.min(), bands.max(), 200) dos, _ = np.histogram(bands.flatten(), bins=energy_bins, density=True) axes[1,0].plot(dos, energy_bins[:-1], 'b-', linewidth=2) axes[1,0].fill_betweenx(energy_bins[:-1], dos, alpha=0.3) axes[1,0].axhline(vbm, color='blue', linestyle='--') axes[1,0].axhline(cbm, color='red', linestyle='--') axes[1,0].set_xlabel('Density of States', fontsize=12) axes[1,0].set_ylabel('Energy (a.u.)', fontsize=12) axes[1,0].set_title('Density of States (DOS)', fontsize=14, fontweight='bold') axes[1,0].grid(True, alpha=0.3) # (4) Summary of results axes[1,1].axis('off') summary = f""" 1D Periodic Crystal Band Structure Lattice Constant: {kp.a:.2f} Potential Depth: {kp.V0:.2f} Plane Waves: {kp.N_pw} Band Gap: {gap:.4f} a.u. {gap * 27.211:.2f} eV VBM: {vbm:.4f} a.u. CBM: {cbm:.4f} a.u. Material Type: {'Insulator' if gap > 0.2 else 'Semiconductor'} """ axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace', verticalalignment='center') plt.tight_layout() plt.savefig('band_structure_1d.png', dpi=150, bbox_inches='tight') plt.show() print(f"\nBand gap: {gap * 27.211:.2f} eV") print(f"Material type: {'Insulator' if gap > 0.2 else 'Semiconductor'}")

class MetalSurfaceAdsorption: """Adsorption energy calculation on a metal surface (simplified model)""" def __init__(self, metal='Pt'): self.metal = metal # Work function of the metal (eV) self.work_functions = {'Pt': 5.65, 'Pd': 5.12, 'Cu': 4.65, 'Au': 5.1} self.phi = self.work_functions.get(metal, 5.0) def lennard_jones(self, r, epsilon, sigma): """Lennard-Jones potential""" return 4 * epsilon * ((sigma/r)**12 - (sigma/r)**6) def adsorption_energy(self, z, molecule='CO'): """Compute the adsorption energy""" # LJ parameters (metal-CO) epsilon = 0.5 # eV sigma = 2.5 # Å # Adsorption energy = LJ + electrostatic interaction E_vdw = self.lennard_jones(z, epsilon, sigma) E_elec = -0.3 * self.phi * np.exp(-z / 2.0) return E_vdw + E_elec def find_optimal_height(self): """Search for the optimal adsorption height""" z_range = np.linspace(1.5, 6.0, 200) E_ads = [self.adsorption_energy(z) for z in z_range] min_idx = np.argmin(E_ads) z_opt = z_range[min_idx] E_opt = E_ads[min_idx] return z_opt, E_opt, z_range, E_ads # Comparison of multiple metals metals = ['Pt', 'Pd', 'Cu', 'Au'] results = {} for metal in metals: surface = MetalSurfaceAdsorption(metal) z_opt, E_opt, z_range, E_ads = surface.find_optimal_height() results[metal] = {'z': z_opt, 'E': E_opt, 'curve': (z_range, E_ads)} print(f"{metal}: z = {z_opt:.2f} Å, E_ads = {E_opt:.3f} eV") # Visualization fig, axes = plt.subplots(2, 2, figsize=(12, 10)) # (1) Adsorption energy curves colors = {'Pt': 'blue', 'Pd': 'green', 'Cu': 'orange', 'Au': 'gold'} for metal, data in results.items(): z_range, E_ads = data['curve'] axes[0,0].plot(z_range, E_ads, color=colors[metal], linewidth=2, label=metal) axes[0,0].plot(data['z'], data['E'], 'o', color=colors[metal], markersize=8) axes[0,0].axhline(0, color='k', linestyle='--', alpha=0.3) axes[0,0].set_xlabel('Distance from Surface (Å)', fontsize=12) axes[0,0].set_ylabel('Adsorption Energy (eV)', fontsize=12) axes[0,0].set_title('CO Adsorption on Metal Surfaces', fontsize=14, fontweight='bold') axes[0,0].legend() axes[0,0].grid(True, alpha=0.3) # (2) Comparison of adsorption energies metal_names = list(results.keys()) E_ads_values = [-results[m]['E'] for m in metal_names] # Displayed as positive values bar_colors = [colors[m] for m in metal_names] axes[0,1].bar(metal_names, E_ads_values, color=bar_colors, alpha=0.7) axes[0,1].set_ylabel('|E_ads| (eV)', fontsize=12) axes[0,1].set_title('Adsorption Energy Comparison', fontsize=14, fontweight='bold') axes[0,1].grid(True, alpha=0.3, axis='y') # Highlight the strongest adsorption max_idx = np.argmax(E_ads_values) axes[0,1].text(max_idx, E_ads_values[max_idx] + 0.05, 'Strongest', ha='center', fontsize=10, fontweight='bold') # (3) Work function vs adsorption energy phi_values = [MetalSurfaceAdsorption(m).phi for m in metal_names] axes[1,0].scatter(phi_values, E_ads_values, s=150, c=bar_colors, alpha=0.7) for i, metal in enumerate(metal_names): axes[1,0].annotate(metal, (phi_values[i], E_ads_values[i]), xytext=(5, 5), textcoords='offset points', fontsize=11, fontweight='bold') # Linear fit z = np.polyfit(phi_values, E_ads_values, 1) p = np.poly1d(z) phi_fit = np.linspace(min(phi_values), max(phi_values), 100) axes[1,0].plot(phi_fit, p(phi_fit), 'r--', alpha=0.5, label=f'y = {z[0]:.2f}x + {z[1]:.2f}') axes[1,0].set_xlabel('Work Function (eV)', fontsize=12) axes[1,0].set_ylabel('|E_ads| (eV)', fontsize=12) axes[1,0].set_title('Work Function vs Adsorption Energy', fontsize=14, fontweight='bold') axes[1,0].legend() axes[1,0].grid(True, alpha=0.3) # (4) Catalytic activity summary axes[1,1].axis('off') best_metal = metal_names[max_idx] summary = f""" CO Adsorption on Metal Surfaces Best Catalyst: {best_metal} E_ads: {-results[best_metal]['E']:.3f} eV Height: {results[best_metal]['z']:.2f} Å Work Func: {MetalSurfaceAdsorption(best_metal).phi:.2f} eV Ranking (by |E_ads|): """ for i, m in enumerate(sorted(metal_names, key=lambda x: -results[x]['E']), 1): summary += f" {i}. {m:3s} {-results[m]['E']:.3f} eV\n" axes[1,1].text(0.1, 0.5, summary, fontsize=11, family='monospace', verticalalignment='center') plt.tight_layout() plt.savefig('catalyst_adsorption.png', dpi=150, bbox_inches='tight') plt.show()