化学反応のシミュレーション：分子衝突から反応速度まで

import numpy as npimport matplotlib.pyplot as plt
np.random.seed(0)
# パラメータN_A, N_B = 100, 100   # A分子・B分子の個数L = 10.0              # 箱の1辺の長さT = 1.0               # 温度（速度スケール）dt = 0.05             # 時間刻み
# 位置はランダムに、速度はガウス分布（→速さがマクスウェル・ボルツマン分布になる）def init(n):    pos = np.random.rand(n, 2) * L    vel = np.random.randn(n, 2) * np.sqrt(T)    return pos, vel
posA, velA = init(N_A)posB, velB = init(N_B)
def step(pos, vel):    pos += vel * dt    # 壁で反射    over = pos < 0    vel[over] *= -1    pos[over] = -pos[over]    over = pos > L    vel[over] *= -1    pos[over] = 2*L - pos[over]    return pos, vel
# 1ステップずつ進めて分布を見るfor _ in range(50):    posA, velA = step(posA, velA)    posB, velB = step(posB, velB)
plt.figure(figsize=(6, 6))plt.scatter(posA[:,0], posA[:,1], c="tab:blue", label="A", s=20)plt.scatter(posB[:,0], posB[:,1], c="tab:red",  label="B", s=20)plt.xlim(0, L); plt.ylim(0, L); plt.legend()plt.title("Step 1: 分子はランダムに動き回る")plt.show()

from scipy.spatial import cKDTree
R_REACT = 0.4   # この距離以下で衝突＝反応
def react(posA, velA, posB, velB):    # 高速な近傍探索    tree = cKDTree(posB)    pairs = tree.query_ball_point(posA, r=R_REACT)    keepA = np.ones(len(posA), dtype=bool)    keepB = np.ones(len(posB), dtype=bool)    for i, near in enumerate(pairs):        if near and keepA[i]:            j = next((k for k in near if keepB[k]), None)            if j is not None:                keepA[i] = False                keepB[j] = False    return posA[keepA], velA[keepA], posB[keepB], velB[keepB]
# シミュレーション本体posA, velA = init(N_A)posB, velB = init(N_B)history = []for t in range(400):    posA, velA = step(posA, velA)    posB, velB = step(posB, velB)    posA, velA, posB, velB = react(posA, velA, posB, velB)    history.append(len(posA))
plt.figure(figsize=(8, 4))plt.plot(history)plt.xlabel("時間ステップ"); plt.ylabel("残ったA分子の数")plt.title("Step 2: 反応が進むにつれてA分子が減る")plt.grid(alpha=0.3); plt.show()

def run(T, steps=400):    posA, velA = np.random.rand(N_A, 2)*L, np.random.randn(N_A, 2)*np.sqrt(T)    posB, velB = np.random.rand(N_B, 2)*L, np.random.randn(N_B, 2)*np.sqrt(T)    hist = []    for _ in range(steps):        posA, velA = step(posA, velA)        posB, velB = step(posB, velB)        posA, velA, posB, velB = react(posA, velA, posB, velB)        hist.append(len(posA))    return np.array(hist)
plt.figure(figsize=(9, 5))for T in [0.5, 1.0, 2.0, 4.0]:    plt.plot(run(T), label=f"T={T}")plt.xlabel("時間ステップ"); plt.ylabel("残ったA分子")plt.title("Step 3: 温度が高いほど反応が速い")plt.legend(); plt.grid(alpha=0.3); plt.show()

Ea = 1.5  # 活性化エネルギー
def react_with_Ea(posA, velA, posB, velB, Ea):    tree = cKDTree(posB)    pairs = tree.query_ball_point(posA, r=R_REACT)    keepA = np.ones(len(posA), dtype=bool)    keepB = np.ones(len(posB), dtype=bool)    for i, near in enumerate(pairs):        if not near or not keepA[i]:            continue        for j in near:            if not keepB[j]:                continue            # 相対運動エネルギー            v_rel = velA[i] - velB[j]            E_kin = 0.5 * np.sum(v_rel**2)            if E_kin >= Ea:                keepA[i] = False                keepB[j] = False                break    return posA[keepA], velA[keepA], posB[keepB], velB[keepB]

from scipy.optimize import curve_fit
def exp_decay(t, N0, k):    return N0 * np.exp(-k * t)
Ts = [0.6, 0.8, 1.0, 1.5, 2.0, 3.0]ks = []for T in Ts:    posA, velA = np.random.rand(N_A,2)*L, np.random.randn(N_A,2)*np.sqrt(T)    posB, velB = np.random.rand(N_B,2)*L, np.random.randn(N_B,2)*np.sqrt(T)    hist = []    for _ in range(800):        posA, velA = step(posA, velA)        posB, velB = step(posB, velB)        posA, velA, posB, velB = react_with_Ea(posA, velA, posB, velB, Ea=1.5)        hist.append(len(posA))    t = np.arange(len(hist))    popt, _ = curve_fit(exp_decay, t, hist, p0=[N_A, 0.01])    ks.append(popt[1])
import numpy as npinv_T = 1.0 / np.array(Ts)ln_k = np.log(ks)
plt.figure(figsize=(8, 5))plt.scatter(inv_T, ln_k, s=80)# 直線フィットslope, intercept = np.polyfit(inv_T, ln_k, 1)xs = np.linspace(inv_T.min(), inv_T.max(), 100)plt.plot(xs, slope*xs + intercept, '--', label=f"傾き = {slope:.2f}")plt.xlabel("1 / T"); plt.ylabel("ln(k)")plt.title("Step 5: アレニウスプロット（直線になる）")plt.legend(); plt.grid(alpha=0.3); plt.show()print(f"推定された活性化エネルギー Ea ≈ {-slope:.2f}（設定値 1.5）")