[Python] MD simulation using velocity verlet in Python

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?

Continue reading...