MD simulation using velocity verlet in Python

Question

I'm trying to implement a simple MD simulation in Python (I'm new to this),I'm using LJ potential and force equations along with Verlet method:

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

Now to make sure it works I'm trying to use this on the simplest system, i.e 2 atoms separated by the distance 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()

But I'm getting incorrect values and the plot shows that the atoms are not bouncing at all, in contra they are drifting away from each other!

I have been trying to find where I went wrong for a very long time now but not progressing in any level! Thought that a second eye might help! Any tips/advice?

Answer 1

Your initial force F is negative. You assign that to P1, making it move down and -F (which is positive) to P2, so making P2 move up. That is contrary to the physics.

You have just assigned F and -F to the wrong particles (F should go to the one with positive r). Switch them round, both where you update x_new, v_new and the particle forces p1_F and p2_F.

With your given potential I think the equilibrium displacement should be 2^(1/6).sigma, or about 3.816 in the units you are using here.

Specifically:

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)

Gives