第5章: 統計力学シミュレーション

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from IPython.display import HTML class Ising2D: def __init__(self, L, T, J=1.0): """ 2次元Isingモデル L: 格子サイズ (L x L) T: 温度 J: 相互作用強度 """ self.L = L self.N = L * L self.T = T self.J = J self.beta = 1.0 / T if T > 0 else np.inf # ランダム初期配置 self.spins = np.random.choice([-1, 1], size=(L, L)) # 統計量 self.energy_history = [] self.magnetization_history = [] def energy(self): """系全体のエネルギー""" E = 0 for i in range(self.L): for j in range(self.L): # 最隣接との相互作用（周期境界条件） S = self.spins[i, j] neighbors = ( self.spins[(i+1) % self.L, j] + self.spins[i, (j+1) % self.L] ) E += -self.J * S * neighbors return E def delta_energy(self, i, j): """スピン (i,j) を反転させたときのエネルギー変化""" S = self.spins[i, j] neighbors = ( self.spins[(i+1) % self.L, j] + self.spins[(i-1) % self.L, j] + self.spins[i, (j+1) % self.L] + self.spins[i, (j-1) % self.L] ) return 2 * self.J * S * neighbors def magnetization(self): """磁化""" return np.abs(np.sum(self.spins)) / self.N def metropolis_step(self): """Metropolis 1ステップ（全スピンを1回ずつ試行）""" for _ in range(self.N): # ランダムにスピンを選択 i, j = np.random.randint(0, self.L, size=2) # エネルギー変化 dE = self.delta_energy(i, j) # Metropolis判定 if dE < 0 or np.random.random() < np.exp(-self.beta * dE): self.spins[i, j] *= -1 # スピン反転 def simulate(self, steps, equilibration=1000): """シミュレーション実行""" # 平衡化 for _ in range(equilibration): self.metropolis_step() # 測定 for step in range(steps): self.metropolis_step() if step % 10 == 0: # 10ステップごとに記録 self.energy_history.append(self.energy()) self.magnetization_history.append(self.magnetization()) def statistics(self): """統計量の計算""" E_mean = np.mean(self.energy_history) E_std = np.std(self.energy_history) M_mean = np.mean(self.magnetization_history) M_std = np.std(self.magnetization_history) # 比熱と磁化率 E_array = np.array(self.energy_history) M_array = np.array(self.magnetization_history) C = self.beta**2 * (np.mean(E_array**2) - np.mean(E_array)**2) / self.N chi = self.beta * self.N * (np.mean(M_array**2) - np.mean(M_array)**2) return { 'E_mean': E_mean / self.N, 'E_std': E_std / self.N, 'M_mean': M_mean, 'M_std': M_std, 'C': C, 'chi': chi } # 異なる温度でシミュレーション L = 20 T_range = np.linspace(1.0, 4.0, 16) J = 1.0 T_c_onsager = 2 * J / np.log(1 + np.sqrt(2)) # Onsager解析解 results = [] for T in T_range: print(f"Temperature T = {T:.2f}...") ising = Ising2D(L, T, J) ising.simulate(steps=5000, equilibration=1000) stats = ising.statistics() stats['T'] = T results.append(stats) # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # エネルギー ax1 = axes[0, 0] T_vals = [r['T'] for r in results] E_vals = [r['E_mean'] for r in results] E_err = [r['E_std'] for r in results] ax1.errorbar(T_vals, E_vals, yerr=E_err, fmt='o-', linewidth=2, markersize=6) ax1.axvline(T_c_onsager, color='r', linestyle='--', linewidth=2, label=f'T_c = {T_c_onsager:.2f}') ax1.set_xlabel('Temperature T') ax1.set_ylabel('Energy per spin ⟨E⟩/N') ax1.set_title('エネルギーの温度依存性') ax1.legend() ax1.grid(True, alpha=0.3) # 磁化 ax2 = axes[0, 1] M_vals = [r['M_mean'] for r in results] M_err = [r['M_std'] for r in results] ax2.errorbar(T_vals, M_vals, yerr=M_err, fmt='o-', linewidth=2, markersize=6, color='blue') ax2.axvline(T_c_onsager, color='r', linestyle='--', linewidth=2, label=f'T_c = {T_c_onsager:.2f}') ax2.set_xlabel('Temperature T') ax2.set_ylabel('Magnetization ⟨|m|⟩') ax2.set_title('磁化の温度依存性') ax2.legend() ax2.grid(True, alpha=0.3) # 比熱 ax3 = axes[1, 0] C_vals = [r['C'] for r in results] ax3.plot(T_vals, C_vals, 'o-', linewidth=2, markersize=6, color='green') ax3.axvline(T_c_onsager, color='r', linestyle='--', linewidth=2, label=f'T_c = {T_c_onsager:.2f}') ax3.set_xlabel('Temperature T') ax3.set_ylabel('Specific heat C') ax3.set_title('比熱（T_c でピーク）') ax3.legend() ax3.grid(True, alpha=0.3) # 磁化率 ax4 = axes[1, 1] chi_vals = [r['chi'] for r in results] ax4.plot(T_vals, chi_vals, 'o-', linewidth=2, markersize=6, color='purple') ax4.axvline(T_c_onsager, color='r', linestyle='--', linewidth=2, label=f'T_c = {T_c_onsager:.2f}') ax4.set_xlabel('Temperature T') ax4.set_ylabel('Susceptibility χ') ax4.set_title('磁化率（T_c でピーク）') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_ising_mc_simulation.png', dpi=300, bbox_inches='tight') plt.show() # 数値結果 print("\n=== 2次元Isingモデル モンテカルロシミュレーション ===\n") print(f"格子サイズ: {L} × {L}") print(f"Onsager臨界温度: T_c = {T_c_onsager:.4f}\n") print("温度依存性:") for r in results: print(f"T = {r['T']:.2f}: E/N = {r['E_mean']:.3f}, M = {r['M_mean']:.3f}, C = {r['C']:.3f}, χ = {r['chi']:.2f}")

import numpy as np import matplotlib.pyplot as plt def autocorrelation(data): """自己相関関数の計算""" n = len(data) mean = np.mean(data) variance = np.var(data) # ゼロ平均化 data_centered = data - mean # 自己相関 autocorr = np.correlate(data_centered, data_centered, mode='full') autocorr = autocorr[n-1:] / (variance * np.arange(n, 0, -1)) return autocorr def correlation_time(autocorr): """相関時間の推定（指数減衰フィット）""" # 最初のゼロ交差または十分小さくなる点を探す for i, val in enumerate(autocorr): if val < np.exp(-1) or val < 0: return i return len(autocorr) // 2 # 2次元Isingモデルでの自己相関 L = 20 temperatures = [1.5, T_c_onsager, 3.0] colors = ['blue', 'red', 'green'] fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 異なる温度での自己相関 ax1 = axes[0, 0] for T, color in zip(temperatures, colors): ising = Ising2D(L, T, J=1.0) ising.simulate(steps=10000, equilibration=1000) # 磁化の自己相関 autocorr = autocorrelation(np.array(ising.magnetization_history)) tau = correlation_time(autocorr) ax1.plot(autocorr[:100], color=color, linewidth=2, label=f'T = {T:.2f}, τ ≈ {tau}') ax1.set_xlabel('Time lag (MC steps)') ax1.set_ylabel('Autocorrelation C(t)') ax1.set_title('磁化の自己相関') ax1.legend() ax1.grid(True, alpha=0.3) ax1.set_yscale('log') # エルゴード性の検証（時間平均 vs アンサンブル平均） ax2 = axes[0, 1] T_test = 2.5 num_runs = 20 time_averages = [] for run in range(num_runs): ising = Ising2D(L, T_test, J=1.0) ising.simulate(steps=5000, equilibration=1000) time_avg = np.mean(ising.magnetization_history) time_averages.append(time_avg) ax2.hist(time_averages, bins=15, color='skyblue', edgecolor='black', alpha=0.7) ax2.axvline(np.mean(time_averages), color='r', linestyle='--', linewidth=2, label=f'Mean = {np.mean(time_averages):.3f}') ax2.set_xlabel('Time-averaged magnetization') ax2.set_ylabel('Frequency') ax2.set_title(f'エルゴード性検証（T = {T_test}, {num_runs} runs）') ax2.legend() ax2.grid(True, alpha=0.3) # 統計誤差の評価 ax3 = axes[1, 0] ising_err = Ising2D(L, T_c_onsager, J=1.0) ising_err.simulate(steps=10000, equilibration=1000) M_data = np.array(ising_err.magnetization_history) autocorr_err = autocorrelation(M_data) tau_err = correlation_time(autocorr_err) # ブロック平均法 block_sizes = np.logspace(0, 3, 20, dtype=int) block_errors = [] for block_size in block_sizes: if block_size > len(M_data) // 2: break n_blocks = len(M_data) // block_size blocks = M_data[:n_blocks * block_size].reshape(n_blocks, block_size) block_means = np.mean(blocks, axis=1) block_error = np.std(block_means) / np.sqrt(n_blocks) block_errors.append(block_error) ax3.semilogx(block_sizes[:len(block_errors)], block_errors, 'o-', linewidth=2, markersize=6) ax3.axhline(block_errors[-1], color='r', linestyle='--', linewidth=2, label=f'Converged error ≈ {block_errors[-1]:.4f}') ax3.set_xlabel('Block size') ax3.set_ylabel('Standard error') ax3.set_title('ブロック平均法による統計誤差評価') ax3.legend() ax3.grid(True, alpha=0.3) # 平衡化プロセス ax4 = axes[1, 1] ising_eq = Ising2D(L, 2.0, J=1.0) # 完全に秩序状態から開始 ising_eq.spins = np.ones((L, L)) equilibration_M = [] for step in range(2000): ising_eq.metropolis_step() if step % 10 == 0: equilibration_M.append(ising_eq.magnetization()) ax4.plot(np.arange(len(equilibration_M)) * 10, equilibration_M, 'b-', linewidth=2) ax4.axhline(np.mean(equilibration_M[-50:]), color='r', linestyle='--', linewidth=2, label='Equilibrium value') ax4.set_xlabel('MC steps') ax4.set_ylabel('Magnetization') ax4.set_title('平衡化プロセス（秩序→無秩序）') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_mc_autocorrelation.png', dpi=300, bbox_inches='tight') plt.show() print("\n=== 自己相関と統計誤差 ===\n") print(f"相関時間 τ ≈ {tau_err} MC steps") print(f"有効サンプル数 N_eff = N / (2τ + 1) = {len(M_data)} / {2*tau_err + 1} ≈ {len(M_data) / (2*tau_err + 1):.0f}") print(f"\n統計誤差（ブロック平均法）: σ ≈ {block_errors[-1]:.4f}")

import numpy as np import matplotlib.pyplot as plt class LennardJonesMD: def __init__(self, N, L, T, dt=0.001, epsilon=1.0, sigma=1.0): """ Lennard-Jones分子動力学 N: 粒子数 L: ボックスサイズ T: 温度 dt: 時間ステップ """ self.N = N self.L = L self.T = T self.dt = dt self.epsilon = epsilon self.sigma = sigma # 初期配置（FCC格子） n = int(np.ceil(N**(1/3))) positions = [] for i in range(n): for j in range(n): for k in range(n): if len(positions) < N: positions.append([i, j, k]) self.positions = np.array(positions, dtype=float) * (L / n) # 初期速度（Maxwell分布） self.velocities = np.random.randn(N, 3) * np.sqrt(T) # 重心運動量をゼロに self.velocities -= np.mean(self.velocities, axis=0) # 履歴 self.energy_history = [] self.temperature_history = [] def apply_pbc(self, r): """周期境界条件""" return r - self.L * np.floor(r / self.L) def lennard_jones(self, r): """Lennard-Jonesポテンシャルと力""" r_norm = np.linalg.norm(r) if r_norm < 0.1 * self.sigma: r_norm = 0.1 * self.sigma # 発散回避 sr6 = (self.sigma / r_norm)**6 potential = 4 * self.epsilon * (sr6**2 - sr6) force_magnitude = 24 * self.epsilon * (2 * sr6**2 - sr6) / r_norm force = force_magnitude * r / r_norm return potential, force def compute_forces(self): """全粒子間の力とポテンシャルエネルギー""" forces = np.zeros((self.N, 3)) potential = 0 for i in range(self.N): for j in range(i+1, self.N): r_ij = self.positions[i] - self.positions[j] r_ij = self.apply_pbc(r_ij) U_ij, F_ij = self.lennard_jones(r_ij) forces[i] += F_ij forces[j] -= F_ij potential += U_ij return forces, potential def velocity_verlet_step(self): """Velocity Verletアルゴリズム""" # 現在の力 forces, potential = self.compute_forces() # 位置更新 self.positions += self.velocities * self.dt + 0.5 * forces * self.dt**2 self.positions = self.apply_pbc(self.positions) # 新しい力 new_forces, new_potential = self.compute_forces() # 速度更新 self.velocities += 0.5 * (forces + new_forces) * self.dt # エネルギーと温度 kinetic = 0.5 * np.sum(self.velocities**2) total_energy = kinetic + new_potential temperature = 2 * kinetic / (3 * self.N) return total_energy, temperature def simulate(self, steps): """シミュレーション実行""" for step in range(steps): energy, temp = self.velocity_verlet_step() if step % 10 == 0: self.energy_history.append(energy) self.temperature_history.append(temp) def radial_distribution(self, bins=100): """動径分布関数 g(r)""" r_max = self.L / 2 hist, bin_edges = np.histogram([], bins=bins, range=(0, r_max)) for i in range(self.N): for j in range(i+1, self.N): r_ij = self.positions[i] - self.positions[j] r_ij = self.apply_pbc(r_ij) r = np.linalg.norm(r_ij) if r < r_max: bin_idx = int(r / r_max * bins) if bin_idx < bins: hist[bin_idx] += 1 # 規格化 r_bins = (bin_edges[:-1] + bin_edges[1:]) / 2 dr = r_bins[1] - r_bins[0] shell_volume = 4 * np.pi * r_bins**2 * dr density = self.N / self.L**3 g_r = hist / (shell_volume * density * self.N) return r_bins, g_r # Lennard-Jones MDシミュレーション N = 64 L = 5.0 T = 1.0 dt = 0.001 steps = 5000 md = LennardJonesMD(N, L, T, dt) md.simulate(steps) # 動径分布関数 r_bins, g_r = md.radial_distribution(bins=50) # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # エネルギー保存 ax1 = axes[0, 0] time = np.arange(len(md.energy_history)) * 10 * dt ax1.plot(time, md.energy_history, 'b-', linewidth=1) ax1.set_xlabel('Time') ax1.set_ylabel('Total Energy') ax1.set_title('エネルギー保存（Verlet法）') ax1.grid(True, alpha=0.3) # 温度揺らぎ ax2 = axes[0, 1] ax2.plot(time, md.temperature_history, 'r-', linewidth=1) ax2.axhline(T, color='k', linestyle='--', linewidth=2, label=f'T_target = {T}') ax2.set_xlabel('Time') ax2.set_ylabel('Temperature') ax2.set_title('温度の時間発展') ax2.legend() ax2.grid(True, alpha=0.3) # 動径分布関数 ax3 = axes[1, 0] ax3.plot(r_bins, g_r, 'g-', linewidth=2) ax3.set_xlabel('r / σ') ax3.set_ylabel('g(r)') ax3.set_title('動径分布関数') ax3.grid(True, alpha=0.3) # 粒子配置（2Dプロジェクション） ax4 = axes[1, 1] ax4.scatter(md.positions[:, 0], md.positions[:, 1], s=100, alpha=0.6, c='blue', edgecolors='black') ax4.set_xlim([0, L]) ax4.set_ylim([0, L]) ax4.set_xlabel('x') ax4.set_ylabel('y') ax4.set_title('粒子配置（XY平面）') ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_lennard_jones_md.png', dpi=300, bbox_inches='tight') plt.show() print("\n=== Lennard-Jones 分子動力学 ===\n") print(f"粒子数: {N}") print(f"ボックスサイズ: {L:.2f}") print(f"時間ステップ: {dt}") print(f"総ステップ数: {steps}") print(f"\n平均エネルギー: {np.mean(md.energy_history):.4f}") print(f"平均温度: {np.mean(md.temperature_history):.4f}") print(f"温度標準偏差: {np.std(md.temperature_history):.4f}")

import numpy as np import matplotlib.pyplot as plt class GCMCAdsorption: def __init__(self, L, mu, T, epsilon_wall=-1.0, sigma=1.0): """ Grand Canonical Monte Carlo 吸着シミュレーション L: ボックスサイズ mu: 化学ポテンシャル T: 温度 epsilon_wall: 壁との相互作用 """ self.L = L self.mu = mu self.T = T self.beta = 1.0 / T self.epsilon_wall = epsilon_wall self.sigma = sigma self.particles = [] # 粒子位置のリスト self.N_history = [] def wall_energy(self, z): """壁との相互作用エネルギー（z方向）""" # 簡易的なLennard-Jones型壁ポテンシャル if z < self.sigma: z = self.sigma E = self.epsilon_wall * ((self.sigma / z)**12 - (self.sigma / z)**6) return E def particle_energy(self, pos): """粒子のエネルギー（壁との相互作用のみ）""" z = pos[2] return self.wall_energy(z) def insert_particle(self): """粒子挿入試行""" # ランダム位置 new_pos = np.random.rand(3) * self.L # エネルギー計算 E_new = self.particle_energy(new_pos) # 受理確率 N = len(self.particles) acceptance = min(1, (self.L**3 / (N + 1)) * np.exp(self.beta * (self.mu - E_new))) if np.random.random() < acceptance: self.particles.append(new_pos) return True return False def delete_particle(self): """粒子削除試行""" if len(self.particles) == 0: return False # ランダムに粒子選択 idx = np.random.randint(len(self.particles)) pos = self.particles[idx] # エネルギー計算 E_old = self.particle_energy(pos) # 受理確率 N = len(self.particles) acceptance = min(1, (N / self.L**3) * np.exp(-self.beta * (self.mu - E_old))) if np.random.random() < acceptance: del self.particles[idx] return True return False def move_particle(self): """粒子移動試行""" if len(self.particles) == 0: return False idx = np.random.randint(len(self.particles)) old_pos = self.particles[idx].copy() # ランダム移動 displacement = (np.random.rand(3) - 0.5) * 0.5 new_pos = old_pos + displacement # 周期境界条件 new_pos = new_pos % self.L # エネルギー差 dE = self.particle_energy(new_pos) - self.particle_energy(old_pos) # Metropolis判定 if dE < 0 or np.random.random() < np.exp(-self.beta * dE): self.particles[idx] = new_pos return True return False def gcmc_step(self): """GCMC 1ステップ""" # ランダムに操作を選択 operation = np.random.choice(['insert', 'delete', 'move'], p=[0.33, 0.33, 0.34]) if operation == 'insert': self.insert_particle() elif operation == 'delete': self.delete_particle() else: self.move_particle() def simulate(self, steps, record_interval=10): """シミュレーション実行""" for step in range(steps): self.gcmc_step() if step % record_interval == 0: self.N_history.append(len(self.particles)) def density_profile(self, bins=50): """密度プロファイル（z方向）""" if len(self.particles) == 0: return np.linspace(0, self.L, bins), np.zeros(bins) z_coords = [p[2] for p in self.particles] hist, bin_edges = np.histogram(z_coords, bins=bins, range=(0, self.L)) z_bins = (bin_edges[:-1] + bin_edges[1:]) / 2 # 密度（粒子数/体積） bin_volume = self.L**2 * (self.L / bins) density = hist / bin_volume return z_bins, density # 異なる化学ポテンシャルでGCMC L = 10.0 T = 1.0 epsilon_wall = -2.0 mu_values = [-5.0, -4.0, -3.0, -2.0] results_gcmc = [] for mu in mu_values: print(f"Chemical potential μ = {mu:.2f}...") gcmc = GCMCAdsorption(L, mu, T, epsilon_wall) gcmc.simulate(steps=10000, record_interval=10) z_bins, density = gcmc.density_profile(bins=50) results_gcmc.append({ 'mu': mu, 'N_avg': np.mean(gcmc.N_history[-100:]), 'N_history': gcmc.N_history, 'z_bins': z_bins, 'density': density }) # 可視化 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # 粒子数の時間発展 ax1 = axes[0, 0] colors = ['blue', 'green', 'orange', 'red'] for r, color in zip(results_gcmc, colors): ax1.plot(r['N_history'], color=color, linewidth=1, alpha=0.7, label=f"μ = {r['mu']:.1f}") ax1.set_xlabel('MC steps') ax1.set_ylabel('Number of particles N') ax1.set_title('粒子数の時間発展') ax1.legend() ax1.grid(True, alpha=0.3) # 吸着等温線 ax2 = axes[0, 1] mu_vals = [r['mu'] for r in results_gcmc] N_avg_vals = [r['N_avg'] for r in results_gcmc] ax2.plot(mu_vals, N_avg_vals, 'o-', linewidth=2, markersize=8, color='purple') ax2.set_xlabel('Chemical potential μ') ax2.set_ylabel('Average number of particles ⟨N⟩') ax2.set_title('吸着等温線（GCMC）') ax2.grid(True, alpha=0.3) # 密度プロファイル ax3 = axes[1, 0] for r, color in zip(results_gcmc, colors): ax3.plot(r['z_bins'], r['density'], color=color, linewidth=2, label=f"μ = {r['mu']:.1f}") ax3.set_xlabel('z (distance from wall)') ax3.set_ylabel('Density ρ(z)') ax3.set_title('密度プロファイル（壁近傍で吸着）') ax3.legend() ax3.grid(True, alpha=0.3) # 粒子数分布 ax4 = axes[1, 1] for r, color in zip(results_gcmc, colors): ax4.hist(r['N_history'][-500:], bins=20, alpha=0.5, color=color, label=f"μ = {r['mu']:.1f}, ⟨N⟩ = {r['N_avg']:.1f}") ax4.set_xlabel('Number of particles N') ax4.set_ylabel('Frequency') ax4.set_title('粒子数分布（揺らぎ）') ax4.legend() ax4.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('stat_mech_gcmc_adsorption.png', dpi=300, bbox_inches='tight') plt.show() print("\n=== GCMC 吸着シミュレーション ===\n") print("結果:") for r in results_gcmc: print(f"μ = {r['mu']:.2f}: ⟨N⟩ = {r['N_avg']:.2f}")