不使用已有的分子动力学模拟库，只用numpy、matplotlib等基础库，怎么用python模拟上述题目？（使用周期边界条件）

使用基础Python库如`numpy`和`matplotlib`来模拟上述问题，特别是涉及Lennard-Jones潜在能和刚性球模型的情况，你需要创建一个简单的自定义算法。以下是大致的步骤，但由于这个问题的复杂性，完整代码将非常长，这里提供一个简化的示例框架：

import numpy as np
import matplotlib.pyplot as plt

# 定义Lennard-Jones函数 (简化版)
def lj_potential(r, sigma=3.405e-10, epsilon=119.8):
    return 4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)

# 周期边界条件
def wrap(positions, box_size):
    return positions % box_size

# 粒子类
class Particle:
    def __init__(self, pos, vel, radius):
        self.pos = wrap(pos, box_size=6)
        self.vel = vel
        self.radius = radius

# 生成初始状态 (16个粒子，随机位置和速度)
box_size = 6
sigma = 3.405e-10
epsilon = 119.8
n_particles = 16
radius = 0.5 * sigma  # 刚性球半径假设
positions = np.random.uniform(0, box_size, size=(n_particles, 2))
velocities = np.zeros_like(positions)

# 主循环 (这里只是一个伪代码，真实模拟会包含更多细节)
for t in range(steps):  # 循环步数
    forces = np.zeros_like(positions)
    for i, particle_i in enumerate(Particles):
        for j, particle_j in enumerate(particles):
            if i != j:
                distance = np.linalg.norm(positions[i] - positions[j])
                force = lj_potential(distance, sigma=sigma, epsilon=epsilon) * (distance - sigma) / distance**7
                forces[i] += force
                forces[j] -= force  # 平衡力
    accelerations = forces / masses  # 根据牛顿第二定律
    velocities += accelerations * dt
    positions += velocities * dt

    # 绘制速度分布 (简化版本，仅显示一维速度)
    speed_distribution = np.linalg.norm(velocities, axis=1)
    plt.hist(speed_distribution, bins='auto', density=True)
    plt.xlabel('Speed')
    plt.ylabel('Probability')

# 使用刚性球模型计算平均自由程
masses = np.full(n_particles, 6.69e-26)  # 假设所有粒子质量相等
viscosity = ...  # 实际上需要查表或实验数据
mean_free_path = ...  # 使用Stokes-Einstein公式计算，这里省略了计算过程

# 输出结果
print(f"Mean free path at {temperature_in_K}K is {mean_free_path}.")

注意，这只是一个简化的例子，实际的模拟将更复杂，包括粒子间的碰撞检测、温度控制（例如使用Verlet积分或其他优化算法）、以及处理周期边界条件。此外，计算动态黏度和其他物理参数通常需要额外的数据或者近似方法。