嘿，我仔细看了你的代码，发现几个核心问题导致粒子没正确应用PBC、不反弹且能量不守恒，咱们一个个拆解解决：

1. 最关键的遗漏：未对新位置应用PBC #

你已经定义了PBC函数，但在模拟循环里计算完new_positions后完全没调用它！粒子跑出盒子后不会被“包裹”回内部，自然看不到反弹效果，这是最直接的原因。

2. LJ势与力的计算逻辑错误 #

你的LJ_VF函数直接用位移矢量计算，但忽略了两个关键点：

• 必须先取位移的标量距离来计算势能和力的大小，再结合位移的正负（方向）得到矢量力；
• 用r == 0判断会有数值精度问题，应该用极小阈值（比如1e-10）避免除以0。

3. 力的方向与维度处理不当 #

原代码里直接加减标量力，没有结合位移的正负判断方向，导致粒子受力方向错误，运动行为异常。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

def LJ_VF(d):
    # d是一维位移矢量（带正负，代表方向）
    r_mag = np.abs(d)
    if r_mag < 1e-10:  # 用小阈值避免除以0
        return 0.0, 0.0  
    # 计算LJ势能（标量）
    V = 4 * epsilon * ((sigma/r_mag)**12 - (sigma/r_mag)**6)
    # 计算力的大小，再结合位移方向得到矢量力
    F_mag = 24 * epsilon * (2 * (sigma**12)/(r_mag**13) - (sigma**6)/(r_mag**7))
    F = F_mag * (d / r_mag)  # 力的方向与位移一致
    return V, F

def position_verlet(x, v, f_old, dt, m=1.0):
    # Verlet位置更新：考虑质量m（F=ma）
    return x + v * dt + 0.5 * (f_old/m) * dt**2

def velocity_verlet(v, f_old, f_new, dt, m=1.0):
    # Verlet速度更新
    return v + 0.5 * (f_old + f_new)/m * dt

def PBC(pos, uc_size):
    return pos - uc_size * np.floor(pos / uc_size)

def MID(r1, r2, uc_size):
    d = r2 - r1
    return d - uc_size * np.round(d / uc_size)

def LJ_N(positions, L):
    N = len(positions)
    total_F = np.zeros(N)  
    total_V = 0.0  
    for i in range(N):
        for j in range(i + 1, N): 
            d = MID(positions[i], positions[j], L)   
            V, F = LJ_VF(d)  
            total_V += V  
            total_F[i] += F  # i受到j的作用力
            total_F[j] -= F  # j受到i的反作用力
    return total_V, total_F

def simulate_N_atoms(N_atoms, dt, L, N_steps):
    # 初始化位置：均匀分布在盒子内
    positions = np.linspace(0, L, N_atoms, endpoint=False)
    # 给初始速度加小扰动（全0的话粒子只会缓慢分开，看不到边界反弹）
    velocities = np.random.normal(0, 0.5, N_atoms)

    V_tot, forces = LJ_N(positions, L)

    pos_list, vel_list = [PBC(positions.copy(), L)], [velocities.copy()]
    V_list, K_list, H_list = [V_tot], [0.5 * np.sum(velocities**2)], [V_tot + 0.5 * np.sum(velocities**2)]

    for _ in range(N_steps - 1):
        # 用定义好的Verlet函数更新位置
        new_positions = position_verlet(positions, velocities, forces, dt)
        # 对新位置应用PBC，确保粒子在盒子内
        new_positions = PBC(new_positions, L)
        # 计算新位置下的力和势能
        new_V_tot, new_forces = LJ_N(new_positions, L)
        # 更新速度
        new_velocities = velocity_verlet(velocities, forces, new_forces, dt)

        positions, velocities, forces = new_positions, new_velocities, new_forces
        pos_list.append(positions.copy())
        vel_list.append(velocities.copy())
        V_list.append(new_V_tot)
        K_list.append(0.5 * np.sum(velocities**2))
        H_list.append(new_V_tot + K_list[-1])

    return np.array(pos_list), np.array(vel_list), np.array(V_list), np.array(K_list), np.array(H_list)

# 常量设置
epsilon = 0.0103       
sigma = 3.4  
dt = 0.01 
N_steps = 1000  # 增加步数，方便观察运动
L = 10.0
N_atoms = 5

# 运行模拟并绘图
positions, velocities, V_list, K_list, H_list = simulate_N_atoms(N_atoms, dt, L, N_steps)
time = np.arange(N_steps) * dt

# 粒子位置随时间变化图
plt.figure(figsize=(12,6))
for i in range(N_atoms):
    plt.plot(time, positions[:, i], label=f"Atom {i+1}")
plt.axhline(y=L, color='black', linestyle="--", label="Box Boundary")
plt.axhline(y=0, color='black', linestyle="--")
plt.xlabel("Time (t)")
plt.ylabel("x Position (Å)")
plt.title("Particle Positions Over Time with PBC")
plt.legend()
plt.grid(True)
plt.show()

# 能量守恒验证图
plt.figure(figsize=(12,6))
plt.plot(time, V_list, label="Potential Energy")
plt.plot(time, K_list, label="Kinetic Energy")
plt.plot(time, H_list, label="Total Energy")
plt.xlabel("Time (t)")
plt.ylabel("Energy")
plt.title("Energy Conservation")
plt.legend()
plt.grid(True)
plt.show()

关键修正说明 #

• 强制应用PBC：计算新位置后立刻调用PBC函数，确保粒子始终在[0, L)范围内，实现“边界反弹”效果；
• 修正LJ力方向：结合位移的正负确定力的方向，避免受力逻辑错误；
• 添加初始速度扰动：原代码初始速度全为0，粒子只会缓慢分开，加小随机速度后能直观看到粒子碰撞边界的行为；
• 规范参数传递：将L作为参数传入LJ_N函数，避免依赖全局变量；
• 增加模拟步数：原来的10步太少，改成1000步能清晰观察粒子运动和能量变化。

运行修正后的代码，你会看到粒子在边界处“反弹”（实际是通过PBC包裹到盒子另一侧），且总能量基本保持守恒，波动在数值精度范围内。

备注：内容来源于stack exchange，提问作者Lana.s