第2章: カノニカル集団と分配関数

import numpy as np import matplotlib.pyplot as plt # 物理定数 k_B = 1.380649e-23 # J/K hbar = 1.054572e-34 # J·s # 調和振動子の分配関数 omega = 2e13 # 角振動数 (rad/s) - 典型的な格子振動 def partition_function_harmonic(omega, T, hbar, k_B): """量子調和振動子の分配関数""" beta = 1 / (k_B * T) x = beta * hbar * omega return np.exp(-x/2) / (1 - np.exp(-x)) def average_energy(omega, T, hbar, k_B): """平均エネルギー""" beta = 1 / (k_B * T) x = beta * hbar * omega return hbar * omega * (0.5 + 1/(np.exp(x) - 1)) def heat_capacity(omega, T, hbar, k_B): """熱容量 C = dE/dT""" beta = 1 / (k_B * T) x = beta * hbar * omega exp_x = np.exp(x) return k_B * x**2 * exp_x / (exp_x - 1)**2 # 温度範囲 T_range = np.logspace(0, 3, 200) # 1 K ~ 1000 K # 各量を計算 Z_vals = [partition_function_harmonic(omega, T, hbar, k_B) for T in T_range] E_vals = [average_energy(omega, T, hbar, k_B) for T in T_range] C_vals = [heat_capacity(omega, T, hbar, k_B) for T in T_range] # 特性温度 theta = hbar * omega / k_B print(f"=== 調和振動子の熱力学 ===") print(f"角振動数 ω = {omega:.2e} rad/s") print(f"特性温度 θ = ℏω/k_B = {theta:.2f} K\n") # 可視化 fig, axes = plt.subplots(1, 3, figsize=(16, 5)) # 分配関数 ax1 = axes[0] ax1.loglog(T_range, Z_vals, 'b-', linewidth=2) ax1.axvline(theta, color='r', linestyle='--', linewidth=1.5, label=f'θ = {theta:.1f} K') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Partition function Z') ax1.set_title('分配関数') ax1.legend() ax1.grid(True, alpha=0.3, which='both') # 平均エネルギー ax2 = axes[1] E_0 = hbar * omega / 2 # 零点エネルギー ax2.semilogx(T_range, np.array(E_vals) / (hbar * omega), 'g-', linewidth=2) ax2.axhline(0.5, color='k', linestyle=':', linewidth=1, label='零点エネルギー') ax2.axvline(theta, color='r', linestyle='--', linewidth=1.5, label=f'θ = {theta:.1f} K') ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Average energy / ℏω') ax2.set_title('平均エネルギー') ax2.legend() ax2.grid(True, alpha=0.3) # 熱容量 ax3 = axes[2] ax3.semilogx(T_range, np.array(C_vals) / k_B, 'r-', linewidth=2) ax3.axhline(1.0, color='k', linestyle=':', linewidth=1, label='古典極限') ax3.axvline(theta, color='r', linestyle='--', linewidth=1.5, label=f'θ = {theta:.1f} K') ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('Heat capacity / k_B') ax3.set_title('熱容量') ax3.legend() ax3.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_harmonic_oscillator.png', dpi=300, bbox_inches='tight') plt.show() # 低温・高温極限 print("低温極限 (T << θ):") T_low = theta / 10 E_low = average_energy(omega, T_low, hbar, k_B) C_low = heat_capacity(omega, T_low, hbar, k_B) print(f" T = {T_low:.2f} K") print(f" ≈ {E_low/(hbar*omega):.4f} ℏω (≈ 零点エネルギー)") print(f" C ≈ {C_low/k_B:.6f} k_B (指数的に減衰)\n") print("高温極限 (T >> θ):") T_high = theta * 10 E_high = average_energy(omega, T_high, hbar, k_B) C_high = heat_capacity(omega, T_high, hbar, k_B) print(f" T = {T_high:.2f} K") print(f" ≈ {E_high/(k_B*T_high):.4f} k_B T (古典的)") print(f" C ≈ {C_high/k_B:.4f} k_B (等分配定理)")

import numpy as np import matplotlib.pyplot as plt # モンテカルロ法で等分配定理を検証 k_B = 1.380649e-23 # J/K m = 4.65e-26 # N₂分子の質量 (kg) def monte_carlo_classical_gas(N_particles, T, k_B, m, n_steps=100000): """古典理想気体のモンテカルロシミュレーション""" # 速度の初期化（Maxwell分布から） sigma_v = np.sqrt(k_B * T / m) v = np.random.normal(0, sigma_v, (N_particles, 3)) # エネルギーのサンプリング kinetic_energies = [] for step in range(n_steps): # 運動エネルギー KE = 0.5 * m * np.sum(v**2, axis=1) kinetic_energies.append(np.mean(KE)) # Metropolisステップ（速度の微小変化） if step < n_steps - 1: dv = np.random.normal(0, 0.1 * sigma_v, (N_particles, 3)) v_new = v + dv # エネルギー変化 KE_old = 0.5 * m * np.sum(v**2) KE_new = 0.5 * m * np.sum(v_new**2) delta_E = KE_new - KE_old # Metropolis判定 if delta_E < 0 or np.random.rand() < np.exp(-delta_E / (k_B * T * N_particles)): v = v_new return np.array(kinetic_energies) # シミュレーション設定 N_particles = 1000 temperatures = [100, 300, 500, 800] # K colors = ['blue', 'green', 'orange', 'red'] fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # 平均運動エネルギー vs 温度 ax1 = axes[0] avg_KE_list = [] for T in temperatures: KE_samples = monte_carlo_classical_gas(N_particles, T, k_B, m, n_steps=50000) avg_KE = np.mean(KE_samples[10000:]) # 平衡化後 avg_KE_list.append(avg_KE) theoretical = [1.5 * k_B * T for T in temperatures] ax1.plot(temperatures, avg_KE_list, 'bo-', markersize=8, linewidth=2, label='Simulation') ax1.plot(temperatures, theoretical, 'r--', linewidth=2, label='Theory: (3/2)k_B T') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Average kinetic energy (J)') ax1.set_title('等分配定理の検証') ax1.legend() ax1.grid(True, alpha=0.3) # 時間発展 ax2 = axes[1] T_demo = 300 KE_samples = monte_carlo_classical_gas(N_particles, T_demo, k_B, m, n_steps=10000) ax2.plot(KE_samples, 'b-', alpha=0.7, linewidth=1) ax2.axhline(1.5 * k_B * T_demo, color='r', linestyle='--', linewidth=2, label=f'Theory: (3/2)k_B T = {1.5*k_B*T_demo:.2e} J') ax2.set_xlabel('Monte Carlo step') ax2.set_ylabel('Average kinetic energy (J)') ax2.set_title(f'時間発展（T = {T_demo} K）') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_equipartition.png', dpi=300, bbox_inches='tight') plt.show() # 結果の比較 print("=== 等分配定理の検証 ===\n") for T, avg_KE in zip(temperatures, avg_KE_list): theoretical_KE = 1.5 * k_B * T error = abs(avg_KE - theoretical_KE) / theoretical_KE * 100 print(f"T = {T} K:") print(f" シミュレーション: = {avg_KE:.6e} J") print(f" 理論値: (3/2)k_B T = {theoretical_KE:.6e} J") print(f" 誤差: {error:.2f}%\n")

import numpy as np import matplotlib.pyplot as plt k_B = 1.380649e-23 # J/K def two_level_partition(epsilon, T, k_B): """2準位系の分配関数""" beta = 1 / (k_B * T) return 1 + np.exp(-beta * epsilon) def two_level_energy(epsilon, T, k_B): """平均エネルギー""" beta = 1 / (k_B * T) return epsilon / (np.exp(beta * epsilon) + 1) def schottky_heat_capacity(epsilon, T, k_B): """Schottky比熱""" x = epsilon / (k_B * T) exp_x = np.exp(x) return k_B * x**2 * exp_x / (exp_x + 1)**2 # エネルギー間隔 epsilon_values = [1e-21, 2e-21, 5e-21] # J colors = ['blue', 'green', 'red'] labels = ['ε = 1×10⁻²¹ J', 'ε = 2×10⁻²¹ J', 'ε = 5×10⁻²¹ J'] T_range = np.linspace(1, 1000, 500) fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 占有確率 ax1 = axes[0, 0] for epsilon, color, label in zip(epsilon_values, colors, labels): p_0 = 1 / (1 + np.exp(-epsilon / (k_B * T_range))) # 基底状態 p_1 = 1 - p_0 # 励起状態 ax1.plot(T_range, p_0, '-', color=color, linewidth=2, label=f'{label} (ground)') ax1.plot(T_range, p_1, '--', color=color, linewidth=2, label=f'{label} (excited)') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Occupation probability') ax1.set_title('準位占有確率') ax1.legend(fontsize=8) ax1.grid(True, alpha=0.3) # 平均エネルギー ax2 = axes[0, 1] for epsilon, color, label in zip(epsilon_values, colors, labels): E_avg = [two_level_energy(epsilon, T, k_B) for T in T_range] ax2.plot(T_range, np.array(E_avg) / epsilon, color=color, linewidth=2, label=label) ax2.axhline(0.5, color='k', linestyle=':', linewidth=1, label='高温極限') ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Average energy / ε') ax2.set_title('平均エネルギー') ax2.legend() ax2.grid(True, alpha=0.3) # Schottky比熱 ax3 = axes[1, 0] for epsilon, color, label in zip(epsilon_values, colors, labels): C_vals = [schottky_heat_capacity(epsilon, T, k_B) for T in T_range] ax3.plot(T_range, np.array(C_vals) / k_B, color=color, linewidth=2, label=label) # ピーク位置 T_peak = epsilon / (2 * k_B) C_peak = schottky_heat_capacity(epsilon, T_peak, k_B) ax3.plot(T_peak, C_peak / k_B, 'o', color=color, markersize=8) ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('Heat capacity / k_B') ax3.set_title('Schottky比熱') ax3.legend() ax3.grid(True, alpha=0.3) # エネルギー間隔依存性（ピーク温度） ax4 = axes[1, 1] epsilon_scan = np.linspace(0.5e-21, 10e-21, 50) T_peak_list = epsilon_scan / (2 * k_B) C_peak_list = [schottky_heat_capacity(eps, eps/(2*k_B), k_B) for eps in epsilon_scan] ax4.plot(epsilon_scan / 1e-21, T_peak_list, 'b-', linewidth=2) ax4_twin = ax4.twinx() ax4_twin.plot(epsilon_scan / 1e-21, np.array(C_peak_list) / k_B, 'r-', linewidth=2) ax4.set_xlabel('Energy gap ε (10⁻²¹ J)') ax4.set_ylabel('Peak temperature (K)', color='b') ax4_twin.set_ylabel('Peak heat capacity / k_B', color='r') ax4.set_title('ピーク特性のε依存性') ax4.tick_params(axis='y', labelcolor='b') ax4_twin.tick_params(axis='y', labelcolor='r') ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_two_level_schottky.png', dpi=300, bbox_inches='tight') plt.show() # 数値結果 print("=== 2準位系とSchottky比熱 ===\n") for epsilon, label in zip(epsilon_values, labels): T_peak = epsilon / (2 * k_B) C_peak = schottky_heat_capacity(epsilon, T_peak, k_B) print(f"{label}:") print(f" 特性温度 θ = ε/k_B = {epsilon/k_B:.2f} K") print(f" ピーク温度 T_peak ≈ {T_peak:.2f} K") print(f" ピーク比熱 C_peak = {C_peak/k_B:.4f} k_B\n")

import numpy as np import matplotlib.pyplot as plt R = 8.314 # J/(K·mol) def einstein_heat_capacity(T, theta_E, R): """Einstein固体のモル熱容量""" x = theta_E / T exp_x = np.exp(x) return 3 * R * x**2 * exp_x / (exp_x - 1)**2 # 異なる物質（Einstein温度が異なる） materials = { 'Pb': 88, # K 'Ag': 215, 'Al': 394, 'C (diamond)': 1860 } colors = ['blue', 'green', 'orange', 'red'] T_range = np.logspace(0, 3, 300) # 1 K ~ 1000 K fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Einstein熱容量 ax1 = axes[0] for (material, theta_E), color in zip(materials.items(), colors): C_V = [einstein_heat_capacity(T, theta_E, R) for T in T_range] ax1.semilogx(T_range, C_V, color=color, linewidth=2, label=f'{material} (θ_E={theta_E}K)') ax1.axhline(3 * R, color='k', linestyle='--', linewidth=1.5, label='Dulong-Petit: 3R') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Molar heat capacity (J/(K·mol))') ax1.set_title('Einstein固体の熱容量') ax1.legend() ax1.grid(True, alpha=0.3) ax1.set_ylim([0, 30]) # 規格化した熱容量（普遍曲線） ax2 = axes[1] for (material, theta_E), color in zip(materials.items(), colors): T_reduced = T_range / theta_E C_V_reduced = [einstein_heat_capacity(T, theta_E, R) / (3*R) for T in T_range] ax2.semilogx(T_reduced, C_V_reduced, color=color, linewidth=2, label=material) ax2.axhline(1.0, color='k', linestyle='--', linewidth=1.5, label='高温極限') ax2.set_xlabel('Reduced temperature T/θ_E') ax2.set_ylabel('C_V / (3R)') ax2.set_title('規格化熱容量（普遍曲線）') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_einstein_solid.png', dpi=300, bbox_inches='tight') plt.show() # 数値結果 print("=== Einstein固体モデル ===\n") print("Dulong-Petitの法則: C_V = 3R = {:.3f} J/(K·mol)\n".format(3*R)) T_test = 300 # K print(f"T = {T_test} K での熱容量:") for material, theta_E in materials.items(): C_V = einstein_heat_capacity(T_test, theta_E, R) print(f" {material:15s}: C_V = {C_V:.3f} J/(K·mol) ({C_V/(3*R)*100:.1f}% of 3R)") print("\n低温での振る舞い（T = 10 K）:") T_low = 10 for material, theta_E in materials.items(): C_V_low = einstein_heat_capacity(T_low, theta_E, R) print(f" {material:15s}: C_V = {C_V_low:.6f} J/(K·mol)")

import numpy as np import matplotlib.pyplot as plt from scipy.integrate import quad R = 8.314 # J/(K·mol) def debye_integral(x): """Debye積分: D(x) = (3/x³) ∫₀ˣ (t⁴ e^t)/(e^t - 1)² dt""" if x < 1e-10: return 0 def integrand(t): if t < 1e-10: return 0 exp_t = np.exp(t) return (t**4 * exp_t) / (exp_t - 1)**2 integral_value, _ = quad(integrand, 0, x, limit=100) return 3 * integral_value / x**3 def debye_heat_capacity(T, theta_D, R): """Debyeモデルの熱容量""" x = theta_D / T D_x = debye_integral(x) return 3 * R * D_x def debye_T3_law(T, theta_D, R): """低温極限のT³則""" return (12 * np.pi**4 / 5) * R * (T / theta_D)**3 # 物質のDebye温度 materials_debye = { 'Pb': 105, # K 'Ag': 225, 'Cu': 343, 'Al': 428, 'C (diamond)': 2230 } colors = ['blue', 'green', 'orange', 'red', 'purple'] T_range = np.logspace(0, 3, 100) fig, axes = plt.subplots(1, 2, figsize=(14, 6)) # Debye熱容量 ax1 = axes[0] for (material, theta_D), color in zip(materials_debye.items(), colors): C_V = [debye_heat_capacity(T, theta_D, R) for T in T_range] ax1.semilogx(T_range, C_V, '-', color=color, linewidth=2, label=f'{material} (θ_D={theta_D}K)') ax1.axhline(3 * R, color='k', linestyle='--', linewidth=1.5, label='Dulong-Petit') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Molar heat capacity (J/(K·mol))') ax1.set_title('Debye固体の熱容量') ax1.legend(fontsize=9) ax1.grid(True, alpha=0.3) # T³則の検証（低温） ax2 = axes[1] T_low = np.linspace(1, 50, 100) material_test = 'Al' theta_D_test = materials_debye[material_test] C_debye = [debye_heat_capacity(T, theta_D_test, R) for T in T_low] C_T3 = [debye_T3_law(T, theta_D_test, R) for T in T_low] ax2.plot(T_low**3, C_debye, 'b-', linewidth=2, label='Debyeモデル') ax2.plot(T_low**3, C_T3, 'r--', linewidth=2, label='T³則（低温極限）') ax2.set_xlabel('T³ (K³)') ax2.set_ylabel('Heat capacity (J/(K·mol))') ax2.set_title(f'T³則の検証（{material_test}, θ_D={theta_D_test}K）') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_debye_model.png', dpi=300, bbox_inches='tight') plt.show() # EinsteinモデルとDebyeモデルの比較 print("=== EinsteinモデルとDebyeモデルの比較 ===\n") print("物質: Al (θ_E ≈ θ_D ≈ 400 K)\n") theta_test = 400 T_values = [10, 50, 100, 300] from scipy.integrate import quad as scipy_quad for T in T_values: C_einstein = einstein_heat_capacity(T, theta_test, R) C_debye = debye_heat_capacity(T, theta_test, R) print(f"T = {T} K:") print(f" Einstein: C_V = {C_einstein:.4f} J/(K·mol)") print(f" Debye: C_V = {C_debye:.4f} J/(K·mol)") print(f" 差: {abs(C_einstein - C_debye):.4f} J/(K·mol)\n") print("低温でのT³則:") T_low_test = 10 C_debye_low = debye_heat_capacity(T_low_test, theta_test, R) C_T3_low = debye_T3_law(T_low_test, theta_test, R) print(f"T = {T_low_test} K:") print(f" Debye完全計算: {C_debye_low:.6f} J/(K·mol)") print(f" T³近似: {C_T3_low:.6f} J/(K·mol)") print(f" 相対誤差: {abs(C_debye_low - C_T3_low)/C_debye_low * 100:.2f}%")

import numpy as np import matplotlib.pyplot as plt # 2原子分子の回転エネルギー # E_J = B J(J+1), J = 0, 1, 2, ... # B = ℏ²/(2I) : 回転定数 # 縮退度: g_J = 2J + 1 k_B = 1.380649e-23 # J/K hbar = 1.054572e-34 # J·s # H₂分子の回転定数 B_H2 = 85.4 # K（B/k_B の単位で） def rotational_energy(J, B): """回転エネルギー準位（温度単位）""" return B * J * (J + 1) def rotational_partition_function(T, B, J_max=50): """回転分配関数""" Z_rot = 0 for J in range(J_max): g_J = 2 * J + 1 E_J = rotational_energy(J, B) Z_rot += g_J * np.exp(-E_J / T) return Z_rot def rotational_heat_capacity(T, B, J_max=50): """回転熱容量""" Z = rotational_partition_function(T, B, J_max) # の計算 E_avg = 0 for J in range(J_max): g_J = 2 * J + 1 E_J = rotational_energy(J, B) E_avg += g_J * E_J * np.exp(-E_J / T) E_avg /= Z # の計算 E2_avg = 0 for J in range(J_max): g_J = 2 * J + 1 E_J = rotational_energy(J, B) E2_avg += g_J * E_J**2 * np.exp(-E_J / T) E2_avg /= Z # C = d/dT = ( - ²) / (k_B T²) C_rot = (E2_avg - E_avg**2) / T**2 return C_rot * k_B # 温度範囲 T_range = np.linspace(10, 1000, 200) # 準位占有数 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 回転準位のエネルギーと縮退度 ax1 = axes[0, 0] J_values = np.arange(0, 20) E_values = [rotational_energy(J, B_H2) for J in J_values] g_values = [2*J + 1 for J in J_values] ax1.stem(J_values, E_values, basefmt=' ', linefmt='b-', markerfmt='bo') for J, E, g in zip(J_values[:10], E_values[:10], g_values[:10]): ax1.text(J + 0.3, E, f'g={g}', fontsize=8) ax1.set_xlabel('Rotational quantum number J') ax1.set_ylabel('Energy E/k_B (K)') ax1.set_title(f'回転準位（H₂, B = {B_H2} K）') ax1.grid(True, alpha=0.3) # 温度での占有分布 ax2 = axes[0, 1] temperatures_demo = [50, 100, 300, 500] colors_demo = ['blue', 'green', 'orange', 'red'] for T, color in zip(temperatures_demo, colors_demo): Z = rotational_partition_function(T, B_H2) populations = [] for J in J_values: g_J = 2 * J + 1 E_J = rotational_energy(J, B_H2) p_J = g_J * np.exp(-E_J / T) / Z populations.append(p_J) ax2.plot(J_values, populations, 'o-', color=color, linewidth=2, markersize=4, label=f'T = {T} K') ax2.set_xlabel('J') ax2.set_ylabel('Population') ax2.set_title('回転準位の占有分布') ax2.legend() ax2.grid(True, alpha=0.3) # 分配関数 ax3 = axes[1, 0] Z_values = [rotational_partition_function(T, B_H2) for T in T_range] ax3.plot(T_range, Z_values, 'b-', linewidth=2) ax3.axvline(B_H2, color='r', linestyle='--', linewidth=1.5, label=f'θ_rot = {B_H2} K') ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('Partition function Z_rot') ax3.set_title('回転分配関数') ax3.legend() ax3.grid(True, alpha=0.3) # 回転熱容量 ax4 = axes[1, 1] C_rot_values = [rotational_heat_capacity(T, B_H2) / k_B for T in T_range] ax4.plot(T_range, C_rot_values, 'g-', linewidth=2) ax4.axhline(1.0, color='k', linestyle='--', linewidth=1.5, label='古典極限: k_B') ax4.axvline(B_H2, color='r', linestyle='--', linewidth=1.5, label=f'θ_rot = {B_H2} K') ax4.set_xlabel('Temperature (K)') ax4.set_ylabel('C_rot / k_B') ax4.set_title('回転熱容量') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_rotation.png', dpi=300, bbox_inches='tight') plt.show() print("=== 2原子分子の回転（H₂）===") print(f"回転定数 B = {B_H2} K\n") for T in [50, 100, 300, 500]: Z_rot = rotational_partition_function(T, B_H2) C_rot = rotational_heat_capacity(T, B_H2) print(f"T = {T} K:") print(f" Z_rot = {Z_rot:.4f}") print(f" C_rot = {C_rot/k_B:.4f} k_B\n") print("高温極限（T >> θ_rot）:") print(" Z_rot ≈ T/θ_rot") print(" C_rot ≈ k_B (古典的)")

import numpy as np import matplotlib.pyplot as plt k_B = 1.380649e-23 # J/K h = 6.62607015e-34 # J·s hbar = h / (2 * np.pi) m_N2 = 4.65e-26 # N₂分子の質量 (kg) # N₂分子のパラメータ B_N2 = 2.88 # K（回転定数） theta_vib_N2 = 3374 # K（振動温度） def translational_partition_per_molecule(V, T, m, h, k_B): """並進分配関数（1分子あたり）""" lambda_T = h / np.sqrt(2 * np.pi * m * k_B * T) # 熱的de Broglie波長 return V / lambda_T**3 def rotational_partition(T, B): """回転分配関数（高温近似）""" return T / B def vibrational_partition(T, theta_vib): """振動分配関数""" return 1 / (1 - np.exp(-theta_vib / T)) def total_partition_per_molecule(V, T, m, B, theta_vib, h, k_B): """全分配関数 = z_trans × z_rot × z_vib""" z_trans = translational_partition_per_molecule(V, T, m, h, k_B) z_rot = rotational_partition(T, B) z_vib = vibrational_partition(T, theta_vib) return z_trans * z_rot * z_vib # 標準状態 T_std = 298.15 # K V_molar = 0.0224 # m³/mol N_A = 6.022e23 V_per_molecule = V_molar / N_A # 温度範囲 T_range = np.linspace(100, 1000, 200) # 各寄与を計算 z_trans_list = [translational_partition_per_molecule(V_per_molecule, T, m_N2, h, k_B) for T in T_range] z_rot_list = [rotational_partition(T, B_N2) for T in T_range] z_vib_list = [vibrational_partition(T, theta_vib_N2) for T in T_range] z_total_list = [total_partition_per_molecule(V_per_molecule, T, m_N2, B_N2, theta_vib_N2, h, k_B) for T in T_range] # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 各寄与の温度依存性 ax1 = axes[0, 0] ax1.semilogy(T_range, z_trans_list, 'b-', linewidth=2, label='Translational') ax1.semilogy(T_range, z_rot_list, 'g-', linewidth=2, label='Rotational') ax1.semilogy(T_range, z_vib_list, 'r-', linewidth=2, label='Vibrational') ax1.set_xlabel('Temperature (K)') ax1.set_ylabel('Partition function (per molecule)') ax1.set_title('各自由度の分配関数') ax1.legend() ax1.grid(True, alpha=0.3, which='both') # 全分配関数 ax2 = axes[0, 1] ax2.semilogy(T_range, z_total_list, 'k-', linewidth=2) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Total partition function') ax2.set_title('全分配関数（z_trans × z_rot × z_vib）') ax2.grid(True, alpha=0.3, which='both') # 寄与の割合 ax3 = axes[1, 0] log_z_trans = np.log(z_trans_list) log_z_rot = np.log(z_rot_list) log_z_vib = np.log(z_vib_list) log_z_total = np.log(z_total_list) fraction_trans = log_z_trans / log_z_total fraction_rot = log_z_rot / log_z_total fraction_vib = log_z_vib / log_z_total ax3.plot(T_range, fraction_trans, 'b-', linewidth=2, label='Translational') ax3.plot(T_range, fraction_rot, 'g-', linewidth=2, label='Rotational') ax3.plot(T_range, fraction_vib, 'r-', linewidth=2, label='Vibrational') ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('Fraction of ln(z)') ax3.set_title('寄与の割合') ax3.legend() ax3.grid(True, alpha=0.3) # 熱容量への寄与 ax4 = axes[1, 1] C_trans = 1.5 * np.ones_like(T_range) # (3/2) k_B C_rot = 1.0 * np.ones_like(T_range) # k_B（高温） C_vib_list = [] for T in T_range: x = theta_vib_N2 / T exp_x = np.exp(x) C_vib = x**2 * exp_x / (exp_x - 1)**2 C_vib_list.append(C_vib) C_total = C_trans + C_rot + np.array(C_vib_list) ax4.plot(T_range, C_trans, 'b--', linewidth=2, label='Translational') ax4.plot(T_range, C_rot, 'g--', linewidth=2, label='Rotational') ax4.plot(T_range, C_vib_list, 'r--', linewidth=2, label='Vibrational') ax4.plot(T_range, C_total, 'k-', linewidth=2, label='Total') ax4.axhline(3.5, color='gray', linestyle=':', linewidth=1.5, label='(7/2)k_B') ax4.set_xlabel('Temperature (K)') ax4.set_ylabel('C_V / k_B') ax4.set_title('熱容量への寄与（N₂）') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_ideal_gas_partition.png', dpi=300, bbox_inches='tight') plt.show() # 数値結果 print("=== N₂理想気体の分配関数（標準状態）===") print(f"T = {T_std} K, V = {V_molar} m³/mol\n") z_trans_std = translational_partition_per_molecule(V_per_molecule, T_std, m_N2, h, k_B) z_rot_std = rotational_partition(T_std, B_N2) z_vib_std = vibrational_partition(T_std, theta_vib_N2) z_total_std = z_trans_std * z_rot_std * z_vib_std print(f"並進: z_trans = {z_trans_std:.6e}") print(f"回転: z_rot = {z_rot_std:.2f}") print(f"振動: z_vib = {z_vib_std:.6f}") print(f"全体: z_total = {z_total_std:.6e}\n") print("熱容量（300 K）:") C_trans_val = 1.5 C_rot_val = 1.0 T_300 = 300 x_300 = theta_vib_N2 / T_300 C_vib_val = x_300**2 * np.exp(x_300) / (np.exp(x_300) - 1)**2 C_total_val = C_trans_val + C_rot_val + C_vib_val print(f" C_trans = {C_trans_val:.2f} k_B") print(f" C_rot = {C_rot_val:.2f} k_B") print(f" C_vib = {C_vib_val:.4f} k_B （振動はほぼ凍結）") print(f" C_total ≈ {C_total_val:.2f} k_B") print(f" 実測値: C_V ≈ 2.5 k_B （並進+回転のみ）")

import numpy as np import matplotlib.pyplot as plt # 非調和ポテンシャル: V(x) = (1/2)kx² - gx³ + fx⁴ # gが非調和性パラメータ k_B = 1.380649e-23 # J/K def anharmonic_potential(x, k, g, f): """非調和ポテンシャル""" return 0.5 * k * x**2 - g * x**3 + f * x**4 def thermal_expansion_coefficient(T, k, g, theta, k_B): """線熱膨張係数（摂動理論）""" # α ∝ g/k × (dE/dT) # 簡易モデル x = theta / T if x > 50: # 極低温 return 0 exp_x = np.exp(x) dE_dT = k_B * x**2 * exp_x / (exp_x - 1)**2 alpha = (g / k) * dE_dT / (k_B * T) return alpha # パラメータ（典型的な固体） k = 100 # N/m（バネ定数） g = 10 # 非調和パラメータ f = 1 # 4次項 theta = 300 # K（Einstein温度） x_range = np.linspace(-0.2, 0.3, 1000) fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # ポテンシャルの形状 ax1 = axes[0, 0] V_harmonic = 0.5 * k * x_range**2 V_anharmonic = anharmonic_potential(x_range, k, g, f) ax1.plot(x_range, V_harmonic, 'b--', linewidth=2, label='調和ポテンシャル') ax1.plot(x_range, V_anharmonic, 'r-', linewidth=2, label='非調和ポテンシャル') ax1.axvline(0, color='k', linestyle=':', linewidth=1) ax1.axhline(0, color='k', linestyle=':', linewidth=1) ax1.set_xlabel('Displacement x') ax1.set_ylabel('Potential V(x)') ax1.set_title('調和・非調和ポテンシャル') ax1.legend() ax1.grid(True, alpha=0.3) ax1.set_ylim([-0.5, 2]) # 平衡位置の温度依存性（定性的） ax2 = axes[0, 1] T_vals = np.linspace(1, 1000, 100) # 非調和性により高温で平衡位置がシフト x_eq_shift = [0.01 * (T / theta) * (g / k) if T > theta/10 else 0 for T in T_vals] ax2.plot(T_vals, x_eq_shift, 'r-', linewidth=2) ax2.set_xlabel('Temperature (K)') ax2.set_ylabel('Equilibrium shift') ax2.set_title('平衡位置の温度依存性（熱膨張）') ax2.grid(True, alpha=0.3) # 熱膨張係数 ax3 = axes[1, 0] alpha_vals = [thermal_expansion_coefficient(T, k, g, theta, k_B) for T in T_vals] ax3.plot(T_vals, alpha_vals, 'b-', linewidth=2) ax3.set_xlabel('Temperature (K)') ax3.set_ylabel('Thermal expansion coefficient (a.u.)') ax3.set_title('線熱膨張係数の温度依存性') ax3.grid(True, alpha=0.3) # Grüneisenパラメータとの関係 ax4 = axes[1, 1] # Grüneisenパラメータ γ = V/C_V × α × B_T # 典型的な値 gamma_typical = 2.0 # 無次元 text_str = f""" Grüneisenパラメータ γ: γ = (V/C_V) × α × B_T 典型的な値: γ ≈ {gamma_typical} 意味: - 格子振動の体積依存性 - 熱膨張と圧縮率の関係 - 非調和性の尺度 """ ax4.text(0.1, 0.5, text_str, fontsize=12, verticalalignment='center', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)) ax4.axis('off') ax4.set_title('Grüneisenパラメータ') plt.tight_layout() plt.savefig('stat_mech_anharmonic_thermal_expansion.png', dpi=300, bbox_inches='tight') plt.show() print("=== 非調和効果と熱膨張 ===\n") print("非調和ポテンシャル: V(x) = (1/2)kx² - gx³ + fx⁴\n") print(f"パラメータ:") print(f" k = {k} N/m") print(f" g = {g} （非調和パラメータ）") print(f" θ = {theta} K\n") print("熱膨張のメカニズム:") print(" 1. 非調和性（g≠0）により、ポテンシャルが非対称") print(" 2. 高温で振動振幅が増大") print(" 3. 非対称性により平衡位置がシフト → 熱膨張\n") T_test = 300 alpha_test = thermal_expansion_coefficient(T_test, k, g, theta, k_B) print(f"T = {T_test} K での熱膨張係数:") print(f" α ≈ {alpha_test:.6e} K⁻¹ (任意単位)\n") print("実際の物質の線熱膨張係数（300 K付近）:") print(" Al: α ≈ 23 × 10⁻⁶ K⁻¹") print(" Cu: α ≈ 17 × 10⁻⁶ K⁻¹") print(" Fe: α ≈ 12 × 10⁻⁶ K⁻¹")