Python中使用velocity verlet进行MD仿真

我试图在Python中实现一个简单的MD模拟(我是新手),我使用LJ势和力方程以及Verlet方法:

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

现在,为了确保它有效,我试图在最简单的系统上使用它,即2个原子相距4:

#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()

但我得到的值不正确,图表显示原子根本没有反弹,相反,它们正在相互漂移!

很长一段时间以来,我一直在努力找出我做错的地方,但在任何层面上都没有进步!我想第二只眼睛可能会有所帮助!有什么建议吗?