🔬 Chapter 2: Quantum Theory of Atoms and Molecules

import numpy as np import matplotlib.pyplot as plt from scipy.special import sph_harm, genlaguerre, factorial # Physical constants (SI units) hbar = 1.054571817e-34 # J·s m_e = 9.1093837015e-31 # kg e = 1.602176634e-19 # C epsilon_0 = 8.8541878128e-12 # F/m a_0 = 4 * np.pi * epsilon_0 * hbar**2 / (m_e * e**2) # Bohr radius # Atomic units (simplified) a_0_au = 1.0 # Bohr Ry = 13.6 # eV (Rydberg constant) def hydrogen_energy(n): """Energy level of the hydrogen atom (eV)""" return -Ry / n**2 def radial_wavefunction(r, n, l, a_0=1.0): """Radial wavefunction R_nl(r) (atomic units)""" rho = 2 * r / (n * a_0) # Laguerre polynomial L = genlaguerre(n - l - 1, 2*l + 1) # Normalization constant norm = np.sqrt((2/(n*a_0))**3 * factorial(n - l - 1) / (2*n*factorial(n + l))) R_nl = norm * np.exp(-rho/2) * rho**l * L(rho) return R_nl def radial_probability(r, n, l, a_0=1.0): """Radial probability density r² |R_nl(r)|²""" R = radial_wavefunction(r, n, l, a_0) return r**2 * R**2 # Radial coordinate r = np.linspace(0, 30, 500) # in units of Bohr # Plot fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Radial wavefunctions ax1 = axes[0, 0] orbitals = [(1, 0, '1s'), (2, 0, '2s'), (2, 1, '2p'), (3, 0, '3s'), (3, 1, '3p'), (3, 2, '3d')] for n, l, label in orbitals[:5]: R = radial_wavefunction(r, n, l, a_0_au) ax1.plot(r, R, linewidth=2, label=label) ax1.set_xlabel('r (Bohr radii)') ax1.set_ylabel('R_nl(r)') ax1.set_title('Radial wavefunctions') ax1.legend() ax1.grid(True, alpha=0.3) ax1.axhline(0, color='k', linewidth=0.5) # Radial probability density ax2 = axes[0, 1] for n, l, label in orbitals[:5]: prob = radial_probability(r, n, l, a_0_au) ax2.plot(r, prob, linewidth=2, label=label) ax2.set_xlabel('r (Bohr radii)') ax2.set_ylabel('r² |R_nl(r)|²') ax2.set_title('Radial probability density') ax2.legend() ax2.grid(True, alpha=0.3) # Energy level diagram ax3 = axes[1, 0] n_max = 6 energies = [hydrogen_energy(n) for n in range(1, n_max+1)] degeneracies = [n**2 for n in range(1, n_max+1)] for n, E, g in zip(range(1, n_max+1), energies, degeneracies): ax3.hlines(E, 0, 1, colors='blue', linewidth=2) ax3.text(1.1, E, f'n={n}, E={E:.2f} eV ({g} states)', fontsize=10, va='center') ax3.set_xlim([-0.1, 3]) ax3.set_ylim([min(energies) * 1.1, 0]) ax3.set_ylabel('Energy (eV)') ax3.set_title('Energy levels of the hydrogen atom') ax3.axhline(0, color='k', linestyle='--', linewidth=1, label='Ionization') ax3.grid(True, alpha=0.3, axis='y') ax3.set_xticks([]) ax3.legend() # Most probable radius ax4 = axes[1, 1] r_max_values = [] n_range = range(1, 10) for n in n_range: for l in range(n): r_fine = np.linspace(0.1, 50, 1000) prob = radial_probability(r_fine, n, l, a_0_au) r_max = r_fine[np.argmax(prob)] r_max_values.append((n, l, r_max)) # Color-code by quantum number colors_l = {0: 'blue', 1: 'red', 2: 'green', 3: 'purple'} markers_l = {0: 'o', 1: 's', 2: '^', 3: 'd'} for l_val in [0, 1, 2]: data = [(n, r_max) for n, l, r_max in r_max_values if l == l_val] if data: ns, r_maxs = zip(*data) label_l = ['s', 'p', 'd', 'f'][l_val] ax4.plot(ns, r_maxs, markers_l[l_val], color=colors_l[l_val], markersize=8, linewidth=2, label=f'{label_l} orbitals') ax4.set_xlabel('Principal quantum number n') ax4.set_ylabel('Most probable radius (Bohr radii)') ax4.set_title('Most probable radius') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('qchem_hydrogen_atom.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("=== Atomic orbitals of the hydrogen atom ===\n") print(f"Bohr radius a₀ = {a_0*1e10:.4f} Å") print(f"Rydberg constant Ry = {Ry} eV\n") print("Energy levels:") for n in range(1, 7): E_n = hydrogen_energy(n) print(f" n = {n}: E = {E_n:.4f} eV ({n**2} degenerate)") print("\nMost probable radius (selected orbitals):") for n, l, r_max in r_max_values[:10]: orbital_name = ['s', 'p', 'd', 'f'][l] print(f" {n}{orbital_name}: r_max = {r_max:.2f} a₀")

import numpy as np import matplotlib.pyplot as plt # Dictionary of atomic numbers and electron configurations electron_configs = { 1: '1s¹', # H 2: '1s²', # He 3: '[He]2s¹', # Li 4: '[He]2s²', # Be 5: '[He]2s²2p¹', # B 6: '[He]2s²2p²', # C 7: '[He]2s²2p³', # N 8: '[He]2s²2p⁴', # O 9: '[He]2s²2p⁵', # F 10: '[He]2s²2p⁶', # Ne 11: '[Ne]3s¹', # Na 18: '[Ne]3s²3p⁶', # Ar 26: '[Ar]3d⁶4s²', # Fe 29: '[Ar]3d¹⁰4s¹', # Cu } element_names = { 1: 'H', 2: 'He', 3: 'Li', 4: 'Be', 5: 'B', 6: 'C', 7: 'N', 8: 'O', 9: 'F', 10: 'Ne', 11: 'Na', 18: 'Ar', 26: 'Fe', 29: 'Cu' } # First ionization energy (eV) ionization_energies = { 1: 13.6, 2: 24.6, 3: 5.4, 4: 9.3, 5: 8.3, 6: 11.3, 7: 14.5, 8: 13.6, 9: 17.4, 10: 21.6, 11: 5.1, 18: 15.8 } def orbital_diagram(config_str): """Orbital diagram for an electron configuration""" # Simple parser (implementation omitted, conceptual) pass # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Periodicity of electron configurations (atomic numbers 1-18) ax1 = axes[0, 0] Z_range = list(range(1, 11)) + [11, 18] configs = [electron_configs.get(Z, '') for Z in Z_range] names = [element_names.get(Z, '') for Z in Z_range] ax1.barh(range(len(Z_range)), Z_range, color='skyblue', edgecolor='black') for i, (Z, name, config) in enumerate(zip(Z_range, names, configs)): ax1.text(Z + 0.5, i, f'{name} ({Z}): {config}', va='center', fontsize=9) ax1.set_yticks(range(len(Z_range))) ax1.set_yticklabels(names) ax1.set_xlabel('Atomic number Z') ax1.set_title('Electron configurations (periodic table order)') ax1.grid(True, alpha=0.3, axis='x') # Periodicity of ionization energy ax2 = axes[0, 1] IE_Z = sorted(ionization_energies.keys()) IE_values = [ionization_energies[Z] for Z in IE_Z] ax2.plot(IE_Z, IE_values, 'o-', linewidth=2, markersize=8, color='red') for Z, IE in zip(IE_Z, IE_values): ax2.text(Z, IE + 0.5, element_names[Z], ha='center', fontsize=9) ax2.set_xlabel('Atomic number Z') ax2.set_ylabel('Ionization energy (eV)') ax2.set_title('First ionization energy') ax2.grid(True, alpha=0.3) # Spin multiplicity (using the C atom as an example) ax3 = axes[1, 0] # C atom: 1s² 2s² 2p² # Two electrons in the 2p orbitals (parallel spins by Hund's rule) orbital_labels = ['2p_x', '2p_y', '2p_z'] spins_up = [1, 1, 0] # number of up spins spins_down = [0, 0, 0] # number of down spins x = np.arange(len(orbital_labels)) width = 0.35 bars_up = ax3.bar(x - width/2, spins_up, width, label='Spin up (↑)', color='blue') bars_down = ax3.bar(x + width/2, spins_down, width, label='Spin down (↓)', color='red') ax3.set_ylabel('Number of electrons') ax3.set_xlabel('2p orbitals') ax3.set_title('2p electron configuration of the C atom (Hund\'s rule)') ax3.set_xticks(x) ax3.set_xticklabels(orbital_labels) ax3.set_ylim([0, 2]) ax3.legend() ax3.grid(True, alpha=0.3, axis='y') # Orbital energy diagram (comparison of hydrogen-like and many-electron atoms) ax4 = axes[1, 1] # Hydrogen-like (n dependence only) E_H = {1: -13.6, 2: -3.4, 3: -1.51} orbitals_H = {'1s': (1, -13.6), '2s': (2, -3.4), '2p': (2, -3.4), '3s': (3, -1.51), '3p': (3, -1.51), '3d': (3, -1.51)} # Many-electron atom (with l dependence, conceptual) E_multi = {'1s': -20, '2s': -6, '2p': -5, '3s': -2.5, '3p': -2, '3d': -0.5} x_pos_H = [0.3, 0.8, 1.0, 1.5, 1.7, 1.9] x_pos_multi = [0.5, 1.0, 1.2, 1.7, 1.9, 2.1] orbital_order = ['1s', '2s', '2p', '3s', '3p', '3d'] for i, orb in enumerate(orbital_order): if orb in orbitals_H: ax4.hlines(orbitals_H[orb][1], x_pos_H[i]-0.1, x_pos_H[i]+0.1, colors='blue', linewidth=2) if orb in E_multi: ax4.hlines(E_multi[orb], x_pos_multi[i]-0.1, x_pos_multi[i]+0.1, colors='red', linewidth=2) ax4.plot([], [], 'b-', linewidth=2, label='Hydrogen-like') ax4.plot([], [], 'r-', linewidth=2, label='Multi-electron') for i, orb in enumerate(orbital_order): ax4.text((x_pos_H[i] + x_pos_multi[i])/2, -22, orb, ha='center', fontsize=10) ax4.set_xlim([0, 2.5]) ax4.set_ylim([-25, 0]) ax4.set_ylabel('Energy (eV)') ax4.set_title('Orbital energy ordering (conceptual diagram)') ax4.legend() ax4.grid(True, alpha=0.3, axis='y') ax4.set_xticks([]) plt.tight_layout() plt.savefig('qchem_electron_config.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("\n=== Electron configurations and the Pauli exclusion principle ===\n") print("Electron configurations of selected elements:") for Z in [1, 2, 6, 7, 8, 10]: print(f" {element_names[Z]} (Z={Z}): {electron_configs[Z]}") print("\nPauli exclusion principle:") print(" - At most 1 electron per quantum state (n, l, m, m_s)") print(" - At most 2 electrons per orbital (a spin pair)") print("\nHund's rule (example: the C atom):") print(" 2p²: ↑ ↑ _ (parallel spins, S=1, ³P ground state)")

import numpy as np import matplotlib.pyplot as plt # Atomic units a_0 = 1.0 # Bohr E_h = 1.0 # Hartree def overlap_integral(R, a_0=1.0): """Overlap integral S(R) (simple approximation)""" S = np.exp(-R/a_0) * (1 + R/a_0 + (R/a_0)**2/3) return S def coulomb_integral(R, a_0=1.0): """Coulomb integral H_AA (simple approximation)""" H_AA = -1/a_0**2 - (1/R) * (1 + 1/R) * np.exp(-2*R/a_0) return H_AA def exchange_integral(R, a_0=1.0): """Exchange integral H_AB (simple approximation)""" S = overlap_integral(R, a_0) H_AB = (-S/R - S/(a_0) * (1 + R/a_0)) * np.exp(-R/a_0) return H_AB def h2plus_energy(R, bonding=True, a_0=1.0): """Energy of H₂⁺ (LCAO approximation)""" S = overlap_integral(R, a_0) H_AA = coulomb_integral(R, a_0) H_AB = exchange_integral(R, a_0) # Nuclear repulsion V_NN = 1 / R if bonding: # Bonding orbital E = (H_AA + H_AB) / (1 + S) + V_NN else: # Antibonding orbital E = (H_AA - H_AB) / (1 - S) + V_NN return E # Internuclear distance R_range = np.linspace(0.5, 10, 200) # Potential energy curves E_bonding = [h2plus_energy(R, bonding=True) for R in R_range] E_antibonding = [h2plus_energy(R, bonding=False) for R in R_range] # Equilibrium internuclear distance and dissociation energy E_bonding_array = np.array(E_bonding) R_eq_idx = np.argmin(E_bonding_array) R_eq = R_range[R_eq_idx] E_eq = E_bonding_array[R_eq_idx] E_dissociation = 0 # Dissociation limit (H + H⁺) D_e = E_dissociation - E_eq # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Potential energy curves ax1 = axes[0, 0] ax1.plot(R_range, E_bonding, 'b-', linewidth=2, label='Bonding (σ_g)') ax1.plot(R_range, E_antibonding, 'r-', linewidth=2, label='Antibonding (σ_u*)') ax1.axhline(0, color='k', linestyle='--', linewidth=1, label='Dissociation limit') ax1.plot(R_eq, E_eq, 'go', markersize=10, label=f'Equilibrium (R={R_eq:.2f} a₀)') ax1.set_xlabel('Internuclear distance R (Bohr)') ax1.set_ylabel('Energy (Hartree)') ax1.set_title('H₂⁺ potential energy curve') ax1.set_ylim([-1.5, 2]) ax1.legend() ax1.grid(True, alpha=0.3) # Shape of the molecular orbitals (1D cross section) ax2 = axes[0, 1] R_vis = 2.0 # internuclear distance (for visualization) x = np.linspace(-5, 5, 200) # Atomic orbitals (1s) phi_A = np.exp(-np.abs(x + R_vis/2)) phi_B = np.exp(-np.abs(x - R_vis/2)) # Molecular orbitals S_vis = overlap_integral(R_vis) psi_bonding = (phi_A + phi_B) / np.sqrt(2 * (1 + S_vis)) psi_antibonding = (phi_A - phi_B) / np.sqrt(2 * (1 - S_vis)) ax2.plot(x, psi_bonding, 'b-', linewidth=2, label='Bonding MO (σ_g)') ax2.plot(x, psi_antibonding, 'r-', linewidth=2, label='Antibonding MO (σ_u*)') ax2.axvline(-R_vis/2, color='k', linestyle=':', linewidth=1, label='Nuclei') ax2.axvline(R_vis/2, color='k', linestyle=':', linewidth=1) ax2.axhline(0, color='k', linewidth=0.5) ax2.set_xlabel('Position (Bohr)') ax2.set_ylabel('Wavefunction ψ(x)') ax2.set_title(f'Molecular orbitals (R = {R_vis} a₀)') ax2.legend() ax2.grid(True, alpha=0.3) # Electron density ax3 = axes[1, 0] rho_bonding = psi_bonding**2 rho_antibonding = psi_antibonding**2 ax3.plot(x, rho_bonding, 'b-', linewidth=2, label='Bonding |ψ|²') ax3.plot(x, rho_antibonding, 'r-', linewidth=2, label='Antibonding |ψ|²') ax3.axvline(-R_vis/2, color='k', linestyle=':', linewidth=1) ax3.axvline(R_vis/2, color='k', linestyle=':', linewidth=1) ax3.set_xlabel('Position (Bohr)') ax3.set_ylabel('Electron density |ψ|²') ax3.set_title('Electron density distribution') ax3.legend() ax3.grid(True, alpha=0.3) # Energy level diagram ax4 = axes[1, 1] # Atomic orbital energy E_1s = -0.5 # 1s energy of the H atom (in Hartree) # Molecular orbital energies (R = R_eq) E_sigma_g = h2plus_energy(R_eq, bonding=True) - 1/R_eq # excluding nuclear repulsion E_sigma_u = h2plus_energy(R_eq, bonding=False) - 1/R_eq # Plot ax4.hlines(E_1s, 0, 0.3, colors='blue', linewidth=3, label='H 1s') ax4.hlines(E_1s, 0.7, 1.0, colors='blue', linewidth=3) ax4.hlines(E_sigma_g, 0.35, 0.65, colors='green', linewidth=3, label='σ_g (bonding)') ax4.hlines(E_sigma_u, 0.35, 0.65, colors='red', linewidth=3, label='σ_u* (antibonding)') # Electron configuration (↑ = 1 electron) ax4.plot(0.5, E_sigma_g, 'o', color='black', markersize=10) ax4.text(0.5, E_sigma_g - 0.15, '↑', fontsize=16, ha='center') ax4.set_xlim([-0.1, 1.1]) ax4.set_ylim([-0.8, 0.2]) ax4.set_ylabel('Energy (Hartree)') ax4.set_title('H₂⁺ molecular orbital energy levels') ax4.set_xticks([0.15, 0.5, 0.85]) ax4.set_xticklabels(['H(A)', 'H₂⁺', 'H(B)']) ax4.legend(loc='upper right') ax4.grid(True, alpha=0.3, axis='y') plt.tight_layout() plt.savefig('qchem_h2plus_molecule.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("\n=== H₂⁺ molecule (LCAO approximation) ===\n") print(f"Equilibrium internuclear distance R_eq = {R_eq:.2f} Bohr = {R_eq * 0.529:.2f} Å") print(f"Equilibrium energy E_eq = {E_eq:.4f} Hartree = {E_eq * 27.2:.2f} eV") print(f"Dissociation energy D_e = {D_e:.4f} Hartree = {D_e * 27.2:.2f} eV") print(f"\nComparison with experimental values:") print(f" R_eq (experiment) ≈ 1.06 Å") print(f" D_e (experiment) ≈ 2.79 eV") print(f"\nThe LCAO approximation is qualitatively correct, but quantitative improvement requires expanding the basis set")

import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize_scalar from scipy.integrate import quad # Hamiltonian of the harmonic oscillator (atomic units) hbar = 1.0 m = 1.0 omega = 1.0 def harmonic_potential(x, omega=1.0): """Harmonic oscillator potential""" return 0.5 * omega**2 * x**2 def trial_wavefunction_gaussian(x, alpha): """Trial function (Gaussian)""" return (alpha / np.pi)**0.25 * np.exp(-alpha * x**2 / 2) def kinetic_energy_expectation(alpha): """Expectation value of kinetic energy ⟨T⟩""" # ⟨T⟩ = ∫ ψ* (-ℏ²/2m d²/dx²) ψ dx = ℏ²α/(4m) return hbar**2 * alpha / (4 * m) def potential_energy_expectation(alpha, omega=1.0): """Expectation value of potential energy ⟨V⟩""" # ⟨V⟩ = ∫ ψ* (1/2 m ω² x²) ψ dx = m ω²/(4α) return m * omega**2 / (4 * alpha) def energy_expectation(alpha, omega=1.0): """Expectation value of total energy E(α)""" T = kinetic_energy_expectation(alpha) V = potential_energy_expectation(alpha, omega) return T + V # Range of the variational parameter alpha_range = np.linspace(0.1, 3.0, 200) E_variational = [energy_expectation(alpha, omega) for alpha in alpha_range] # Optimal parameter result = minimize_scalar(lambda a: energy_expectation(a, omega), bounds=(0.1, 10), method='bounded') alpha_opt = result.x E_opt = result.fun # Exact solution E_exact = 0.5 * hbar * omega # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Energy expectation value vs variational parameter ax1 = axes[0, 0] ax1.plot(alpha_range, E_variational, 'b-', linewidth=2, label='E(α)') ax1.axhline(E_exact, color='r', linestyle='--', linewidth=2, label=f'Exact E₀ = {E_exact:.3f}') ax1.plot(alpha_opt, E_opt, 'go', markersize=10, label=f'Optimum α = {alpha_opt:.3f}') ax1.set_xlabel('Variational parameter α') ax1.set_ylabel('Energy E(α)') ax1.set_title('Variational energy') ax1.legend() ax1.grid(True, alpha=0.3) # Comparison of the trial function and the exact solution ax2 = axes[0, 1] x = np.linspace(-4, 4, 500) psi_trial = trial_wavefunction_gaussian(x, alpha_opt) psi_exact = (omega / np.pi)**0.25 * np.exp(-omega * x**2 / 2) ax2.plot(x, psi_trial, 'b-', linewidth=2, label=f'Trial (α={alpha_opt:.3f})') ax2.plot(x, psi_exact, 'r--', linewidth=2, label='Exact') ax2.set_xlabel('Position x') ax2.set_ylabel('Wavefunction ψ(x)') ax2.set_title('Comparison of wavefunctions') ax2.legend() ax2.grid(True, alpha=0.3) # Kinetic energy and potential energy ax3 = axes[1, 0] T_values = [kinetic_energy_expectation(alpha) for alpha in alpha_range] V_values = [potential_energy_expectation(alpha, omega) for alpha in alpha_range] ax3.plot(alpha_range, T_values, 'b-', linewidth=2, label='⟨T⟩') ax3.plot(alpha_range, V_values, 'r-', linewidth=2, label='⟨V⟩') ax3.plot(alpha_range, E_variational, 'g-', linewidth=2, label='⟨E⟩ = ⟨T⟩ + ⟨V⟩') ax3.axvline(alpha_opt, color='k', linestyle=':', linewidth=1, label=f'α_opt') ax3.set_xlabel('Variational parameter α') ax3.set_ylabel('Energy') ax3.set_title('Energy components') ax3.legend() ax3.grid(True, alpha=0.3) # Verification of the virial theorem ax4 = axes[1, 1] # For the harmonic oscillator, ⟨T⟩ = ⟨V⟩ (virial theorem) T_at_opt = kinetic_energy_expectation(alpha_opt) V_at_opt = potential_energy_expectation(alpha_opt, omega) ratio = np.array(T_values) / np.array(V_values) ax4.plot(alpha_range, ratio, 'purple', linewidth=2, label='⟨T⟩ / ⟨V⟩') ax4.axhline(1, color='r', linestyle='--', linewidth=2, label='Virial theorem (=1)') ax4.axvline(alpha_opt, color='k', linestyle=':', linewidth=1) ax4.set_xlabel('Variational parameter α') ax4.set_ylabel('⟨T⟩ / ⟨V⟩') ax4.set_title('Virial theorem') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('qchem_variational_method.png', dpi=300, bbox_inches='tight') plt.show() # Numerical results print("\n=== Ground-state calculation by the variational method ===\n") print(f"Trial function: ψ(x; α) = (α/π)^(1/4) exp(-αx²/2)") print(f"\nOptimization results:") print(f" Optimal parameter α_opt = {alpha_opt:.6f}") print(f" Variational energy E(α_opt) = {E_opt:.6f}") print(f" Exact energy E_exact = {E_exact:.6f}") print(f" Relative error = {abs(E_opt - E_exact)/E_exact * 100:.6f} %") print(f"\nVerification of the virial theorem:") print(f" ⟨T⟩ = {T_at_opt:.6f}") print(f" ⟨V⟩ = {V_at_opt:.6f}") print(f" ⟨T⟩ / ⟨V⟩ = {T_at_opt / V_at_opt:.6f} ≈ 1")