用速度Verlet算法实现等温过程的模拟。编写Python程序实现6×6=36个粒子在二维平面的运动；初始时刻粒子均匀分布，初始速度随机分布（高斯分布）；粒子间相互作用以L-J势能描述；周期性边界条件；统计粒子的均方根位移，并计算自扩散系数。生成动图

速度Verlet算法，也称为半时间步长法，是一种经典的分子动力学模拟方法，用于精确模拟系统的位置和动量。要实现等温过程的模拟，我们需要维护系统的温度并控制其能量。以下是Python程序的一个简化版本，它实现了6x6个粒子在二维平面上的运动：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 参数设置
N = 36
box_size = 10
dt = 0.01
T = 300  # 温度 (K)
k_B = 1.38e-23  # 开尔文-玻尔兹曼常数
gamma = 5  # 摄尔修斯冷却速率
steps = 100000  # 总步骤数

def LJ_potential(r):
    return 4 * ((sigma / r)**12 - (sigma / r)**6)

def force_on_particle(i, pos, box_size, masses):
    forces = [np.zeros(2) for _ in range(N)]
    for j in range(N):
        if i != j:
            rij = pos[j] - pos[i]
            rij %= box_size
            r = np.linalg.norm(rij)
            f = (-LJ_potential(rij) / r**13 + LJ_potential(rij) / r**7) * masses[j]
            forces[i] += f
    return forces[i]

# 初始化
pos = np.random.uniform(low=-box_size/2, high=box_size/2, size=(N, 2))
vels = norm.rvs(size=(N, 2), scale=np.sqrt(T * k_B / masses))  # 高斯分布
masses = np.ones(N)  # 假设所有粒子质量相等
kBT = T * k_B

# Verlet算法核心
for step in range(steps):
    vel_new = vels + dt * (force_on_particle(np.arange(N), pos, box_size, masses) / masses) + gamma * dt * vels / N
    pos_next = pos + vel_new * dt + 0.5 * dt ** 2 * force_on_particle(np.arange(N), pos + vel_new * dt, box_size, masses)
    
    # 碰撞处理和更新状态
    pos_next = np.mod(pos_next, box_size)  # 周期性边界条件
    energies = LJ_potential(np.abs(pos_next[:, None, :] - pos_next[None, :, :]))  # 计算所有对之间的势能
    kinetic_energy = 0.5 * np.sum(masses * vels ** 2)
    total_energy = kinetic_energy + np.sum(energies)
    factor = np.exp(-dt * (total_energy - kBT * N) / k_BT)
    pos, vels = pos_next * factor, vel_new * factor
    
# 统计和绘图
mean_squared_displacement = np.mean((pos[-1] - pos[0])**2, axis=0)
diffusion_coefficient = 2 * mean_squared_displacement / (steps * dt)

print("Mean Squared Displacement:", mean_squared_displacement)
print("Diffusion Coefficient:", diffusion_coefficient)

plt.scatter(pos[:, 0], pos[:, 1])
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Particle Distribution at the End of Simulation")
plt.show()

这个程序会生成一个粒子分布的图像以及最后一步的平均平方位移和扩散系数。注意，这只是一个基本的示例，实际应用中可能需要更复杂的碰撞检测、热浴或者其他调整来维持恒定的温度。