我正在尝试在Python中实现一个简单的MD模拟（我是新手），我将使用LJ电位和力方程与Verlet方法一起使用。 我昨天发布了另一个有关同一代码的问题，终于感谢您取得了进展！ 这是我的代码：

Question

现在，当我绘制速度时，我可以看到它们在增加，因此在应该的情况下，能量并不能保存！

我一直在试图找到很长时间以来我出错的地方，但在任何层面上都没有进步！认为第二只眼可能会有所帮助！有任何提示/建议吗？

（您的帖子以错误的方式恢复为旧代码，因此不会产生您的情节。 您正在慢慢丢失准确性，因为您使用f_new来计算速度，但是F_New是最后一次计算的位置x_old。基本上，f_new不能反映最新位置。 您应该（按顺序）：

计算新位置（使用当前位置，当前速度，旧力）;

计算新力量（取决于新的位置）；

Answer 1

计算新的速度（取决于当前速度和旧力和新力的平均值）。

以下图由以下代码创建。 （原则上，请确保您包括一个，完整的，运行的代码，包括进口。）

import numpy as np
import matplotlib.pyplot as plt

def LJ_VF(r):
    #r = distance in A
    #Returns V in (eV) and F in (eV/A)
    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 position_verlet(x, v, f_old):
    return x + v * dt + 0.5 * f_old * dt**2

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


#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):
        x1_new = position_verlet(p1_x[-1], p1_v[-1], p1_F[-1] )      # update POSITION
        x2_new = position_verlet(p2_x[-1], p2_v[-1], p2_F[-1] )

        r_new = abs(x2_new - x1_new)
        V_new, F_new = LJ_VF(r_new)                                  # update FORCE

        v1_new = velocity_verlet(p1_v[-1], p1_F[-1], -F_new)         # update VELOCITY
        v2_new = velocity_verlet(p2_v[-1], p2_F[-1],  F_new)

        p1_x.append(x1_new)
        p2_x.append(x2_new)
        p1_v.append(v1_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_v, label="Particle 1", color="red")
plt.plot(time, p2_v, label="Particle 2", color="orange")
plt.xlabel("Time (t)")
plt.ylabel("v_x(t)")
plt.title("Velocities v_x(t)")
plt.legend()
plt.show()