Verlet alghoritm and Lennard Jones simulation in Python

[iwko]

Hello,
I created script which draw plot of Lennard Jones potential for two particles. Now i want to create 2D simulation for 10 particles with verlet alghoritm but I have no idea how to create that. My problem is implementation of verlet alghoritm and calculating forces. I implemented boundaries and drawing of particles but I do not know how to calculate force properly using verlet alghoritm so I can't implement moving properly.Someone have solution of that problem? Any advices would be useful.
Below I present my code for plot.

import numpy as np
import matplotlib.pyplot as plt

sigma = 1
e = 1
r = np.zeros([1001])
v = np.zeros([1001])
t = np.zeros([1001])
a = np.zeros([1001])
dt = 0.01
r[0] = 1

def LJ(r):
    return -24*e*((2/r*(sigma/r)**12)-1/r*(sigma/r)**6)

def verlet():
    for i in range(0,1000):
        a[i] = -LJ(r[i])
        r[i+1] = r[i] + dt*v[i]+0.5*dt**2*a[i]
        a[i+1] = -LJ(r[i+1])
        v[i+1] = v[i] + 0.5*dt*(a[i]+a[i+1])
        t[i+1] = t[i]+dt

verlet()
plt.plot(r,a)
plt.axis([1, 2, -2.5, 2.5 ])
plt.show()

 

[DrClaude]

Start by writing down the equations for the Verlet algorithm.

 

[Arsenic&Lace]

Just calculate the LJ force (F_LJ=-grad(LJ)) and use that to obtain the acceleration term.

EDIT: Moreover, the LJ potential is defined based upon the distance between two particles. The way this is currently written, you don't really have a two body problem, unless your coordinate system places one of the particles at the center, which isn't going to scale with the 10 particles. For 10 particles you might want to define periodic boundary conditions with a fixed box size.