Molecular dynamics: Lennard-Jones forces in 2D system

As I am new to all these stuff I would like to verify whether my procedure is correct.

I need to calculate the Lennard-Jones potential, potential energy of the system and the forces.

The procedure is as follows:

1) Read the input data: side lenght of the square -a, number of particles -N, the parameters $\sigma$, $\epsilon$

2) For $i=1,n$ insert randomly n points in x and y axis.

• 2a) For $j=i+1,n+1$ calculate the difference to square between particles in x and y axis.

• 2b) If $r_{ij} \leq 2^{1/6}\sigma$ go back to 2a) step

• 2c) End of the both loops

3) For $i=1,n$

• 3a) For $j=i+1,+1$ calculate the lennard jones potential $$u_{ij}=4\epsilon[(\sigma/r_{ij})^{12}-(\sigma/r_{ij})^6]$$ where $r_{ij}$ is being calculated in 2a) step
• 3b) calculate the potential energy of the system $$U_N=\sum_{i=1}^N\sum_{j=i+1}^N u_{ij}(r_{ij})$$
• 3c) Calculate the force $$F=\sum_{j=1,j\neq i}^NF_{ij}$$ where $$F_{ij}=(-\nabla_ju_{ij}(r_{ij}))$$
• 3d) End of the both loops

I would really grateful for checking, whether it is even correct. I would also appreciate any books connected to this problem - I am currently using Allend Tildesley, but it is rather difficult for me at start.

• Why do you need step 2b? As far as I can tell, it determines whether two particles i and j are too close to each other, this will lead to high values of the pairwise potential $u_{ij}$ that you cannot neglect. – Juan M. Bello-Rivas Sep 6 '15 at 13:51
• Because the particles cannot overlap. – Beginner in fort Sep 6 '15 at 19:16
• I'm not sure what you mean by "loss" but you don't have to do anything special with those particles. Another thing that is not correct is the loop for $j$ = $i+1$ to $N+1$ (it should be $N$, not $N+1$). Regarding 3b and 3c, keep in mind that you have to compute $u_{ij}$ and $-\nabla u_{ij}$ and accumulate those results into $U_N$ and $F$. The rest of what you wrote looks OK. Good luck! – Juan M. Bello-Rivas Sep 6 '15 at 21:10
• Oh, and don't forget that that $F_{ij} = -F_{ji}$. – Juan M. Bello-Rivas Sep 6 '15 at 21:18
• I misread 2b and thought it was during the simulation instead of being the set up, you are free to regenerate those points of course. Now I think you have everything in place. – Juan M. Bello-Rivas Sep 6 '15 at 22:55