基于Velocity Verlet算法的Python分子动力学（MD）双原子模拟异常问题求助

我是分子动力学模拟的新手，最近尝试用Python实现一个基于Lennard-Jones（LJ）势和Velocity Verlet算法的简单MD模拟，先贴出我写的核心函数代码：

def LJ_VF(r):
    #r = distance in Å
    #Returns V in (eV) and F in (eV/Å)
    V = 4 * epsilon * ( (sigma/r)**(12) - (sigma/r)**6 )
    F = 24 * epsilon * ( 2 * ((sigma**12)/(r**(13))) - ( (sigma**6)/(r**7) )) 
    return V , F

def velocity_verlet(x, v, f_old, f_new):                   #setting m=1 so that a=f
    x_new = x + v * dt + 0.5 * f_old * dt**2  
    v_new = v + 0.5 * (f_old + f_new) * dt  
    return x_new, v_new

为了验证代码正确性，我选用了最简单的双原子系统——两个初始间距为4Å的原子，初始速度都设为0。相关的完整模拟代码如下：

#Constants

epsilon = 0.0103     
sigma = 3.4  
m = 1.0
t0 = 0.0
v0 = 0.0
dt = 0.1 
N = 1000

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 = abs(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):
        r_new = abs(p2_x[-1] - p1_x[-1])  
        V_new, F_new = LJ_VF(r_new)  
        x1_new, v1_new = velocity_verlet(p1_x[-1], p1_v[-1], p1_F[-1], F_new)
        x2_new, v2_new = velocity_verlet(p2_x[-1], p2_v[-1], p2_F[-1], -F_new)

        p1_x.append(x1_new)
        p1_v.append(v1_new)
        p2_x.append(x2_new)
        p2_v.append(v2_new)
        
        p1_F.append(F_new)
        p2_F.append(-F_new)
        
        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)

#Initial conditions

p1_x0 = 0.0
p1_v0 = 0.0  
p2_x0 = 4.0  
p2_v0 = 0.0 

#Plot 

p1_x, p1_v, p2_x, p2_v = simulate_two_atoms(p1_x0, p1_v0, p2_x0, p2_v0)
time = np.arange(N)  

plt.plot(time, p1_x, label="Particle 1", color="blue")
plt.plot(time, p2_x, label="Particle 2", color="green")
plt.xlabel("Time (t)")
plt.ylabel("x (Å)")
plt.title("Particle Positions Over Time (Bouncing Test)")
plt.legend()
plt.grid(True)
plt.show()

但运行后结果完全不符合预期：两个原子没有在LJ势的作用下产生振动或反弹，反而彼此越飘越远！得到的位置-时间曲线显示两条平行的直线，原子间距随时间持续增大，完全没有出现预期的相互作用行为。

我已经反复排查了很久，还是找不到代码中的问题，希望有经验的朋友能帮我看看哪里出错了？

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