基于Velocity Verlet算法的Python双原子LJ势分子动力学模拟能量不守恒问题排查

嗨，看起来你已经在MD模拟的路上迈出了不错的第一步！不过双原子模拟里能量不守恒、速度一直上涨的问题，我帮你揪出两个关键的小错误：

1. 力的方向丢失（绝对值导致符号错误） #

你的LJ_VF函数里用了绝对值计算距离，但力是矢量，它的方向直接决定了原子是被吸引还是排斥。当你取绝对值后，返回的F永远是正的，完全丢失了方向信息：

• 当原子间距r > sigma时，LJ势是吸引力，力的方向应该让原子靠近；
• 当r < sigma时，是排斥力，方向应该让原子远离。

错误的力方向会直接导致能量不守恒，甚至让原子一直加速远离/靠近，完全不符合物理规律。

修复方法：
修改LJ_VF函数，保留相对距离的符号，让返回的力自带正确方向：

def LJ_VF(r):
    # r = 带符号的相对距离（比如p2_x - p1_x），单位Å
    # 返回V（标量，eV）和F（矢量，eV/Å，方向沿r轴）
    r_abs = abs(r)
    V = 4 * epsilon * ( (sigma/r_abs)**(12) - (sigma/r_abs)**6 )
    # 先算力的大小，再乘以单位矢量保留方向
    F_magnitude = 24 * epsilon * ( 2 * ((sigma**12)/(r_abs**(13))) - ( (sigma**6)/(r_abs**(7)) ))
    F = F_magnitude * (r / r_abs)  # 单位矢量赋予正确方向
    return V, F

同时在模拟函数里，不要对相对距离取绝对值：

# 初始化时
r = p2_x0 - p1_x0  # 保留符号的相对距离
V, F = LJ_VF(r)
p1_F.append(F)  # p1受到的力指向p2
p2_F.append(-F) # p2受到反作用力，指向p1

# 循环内也一样
r_new = p2_x[-1] - p1_x[-1]  # 不做绝对值处理
V_new, F_new = LJ_VF(r_new)

2. Velocity Verlet的步骤顺序完全搞反了 #

Velocity Verlet的核心是严格的执行顺序，这是保证能量守恒的关键：

1. 用旧位置的力f_old更新原子到新位置；
2. 用新位置计算新力f_new；
3. 用f_old和f_new共同更新速度。

但你的循环里，先基于旧位置计算了F_new，然后直接用它去更新位置和速度——相当于用旧力当新力来计算，积分误差会被不断放大，最终导致能量完全不守恒。

修复后的模拟函数循环部分：

def simulate_two_atoms(p1_x0, p1_v0, p2_x0, p2_v0):
    p1_x, p2_x = [p1_x0], [p2_x0]
    p1_v, p2_v = [p1_v0], [p2_v0]
    p1_F, p1_V, p2_F, p2_V = [], [], [], []

    # 初始化：用初始位置计算初始力
    r = p2_x0 - p1_x0
    V, F = LJ_VF(r)
    p1_F.append(F)
    p1_V.append(V)
    p2_F.append(-F)
    p2_V.append(V)

    for i in range(N - 1):
        # 第一步：用旧力更新位置
        x1_new = p1_x[-1] + p1_v[-1] * dt + 0.5 * p1_F[-1] * dt**2
        x2_new = p2_x[-1] + p2_v[-1] * dt + 0.5 * p2_F[-1] * dt**2
        
        # 第二步：用新位置计算新力
        r_new = x2_new - x1_new
        V_new, F_new = LJ_VF(r_new)
        new_p1_F = F_new
        new_p2_F = -F_new
        
        # 第三步：用旧力+新力更新速度
        v1_new = p1_v[-1] + 0.5 * (p1_F[-1] + new_p1_F) * dt
        v2_new = p2_v[-1] + 0.5 * (p2_F[-1] + new_p2_F) * dt

        # 保存数据
        p1_x.append(x1_new)
        p2_x.append(x2_new)
        p1_v.append(v1_new)
        p2_v.append(v2_new)
        p1_F.append(new_p1_F)
        p2_F.append(new_p2_F)
        p1_V.append(V_new)
        p2_V.append(V_new)
    
    return np.array(p1_x), np.array(p1_v), np.array(p2_x), np.array(p2_v)

你也可以把Velocity Verlet拆成两个独立函数，让逻辑更清晰：

# 拆分位置和速度更新函数，避免混淆
def update_position(x, v, f_old):
    return x + v * dt + 0.5 * f_old * dt**2

def update_velocity(v, f_old, f_new):
    return v + 0.5 * (f_old + f_new) * dt

修复后的预期效果 #

修正这两个问题后，你再运行模拟，会看到两个原子在平衡位置（~3.4Å）附近做简谐振荡，速度会在正负值之间交替变化，总能量（动能+势能）基本保持恒定——这就符合双原子系统的物理规律啦！

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