```
import pygame
import random
import math
import numpy as np
import matplotlib.pyplot as plt
fx = np.zeros([11])
fy = np.zeros([11])
x = np.zeros([11])
y = np.zeros([11])
x[0] = 11
y[0] = 1
sigma = 1.01
e = 1.1
dt = 0.1
def LJ(x,y):
for i in range(1,10):
for j in range(1,10):
rx = (x[i]-x[j])
ry = (y[i]-y[j])
fx[i] = 24*e*(((2/rx)*sigma/rx**12)-((1/rx)*sigma/rx**6))
fy[i] = 24*e*(((2/ry)*sigma/ry**12)-((1/ry)*sigma/ry**6))
print fx, fy
```

Why am I still getting errors
`RuntimeWarning: divide by zero encountered in double_scalars`

and

```
RuntimeWarning: invalid value encountered in double_scalars
```

The result I am getting is

```
[ 0. nan nan nan nan nan nan nan nan nan 0.] [ 0. nan nan nan nan nan nan nan nan nan 0.]
```

I tried to modify starting x and y but that gave no effect.

### Answer1:

In this code:

```
def LJ(x,y):
for i in range(1,10):
for j in range(1,10):
...
```

If `i == j`

, you are comparing a particle to itself. Try skipping that iteration of the for-loops like so:

```
def LJ(x,y):
for i in range(1,10):
for j in range(1,10):
if i == j:
continue
rx = (x[i]-x[j])
ry = (y[i]-y[j])
fx[i] = 24*e*(((2/rx)*sigma/rx**12)-((1/rx)*sigma/rx**6))
fy[i] = 24*e*(((2/ry)*sigma/ry**12)-((1/ry)*sigma/ry**6))
```

Also, you will need to put in actual values for the x and y lists, as they are all currently 0's. Two particles at exactly the same location exert an infinite amount of force according to that equation, so the divide by 0 is accurate in that scenario.

### Answer2:

You should stay with the physics. The potential function is, reconstructing from your force formula,

```
LJ(r) = 4*e*( (sigma/r)**12 - (sigma/r)**6 )
```

The force corresponding to the potential is, for any potential,

```
F = - LJ'(r)*(rx/r, ry/r)
```

Thus the procedure to compute the combined interaction forces should look like

```
def LJ(x,y):
for i in range(x.size):
for j in range(x.size):
rx = (x[i]-x[j])
ry = (y[i]-y[j])
r = sqrt(rx*rx+ry*ry+1e-12)
f = 24*e/r*( 2*(sigma/r)**12 - (sigma/r)**6 )
fx[i] += f*(rx/r)
fy[i] += f*(ry/r)
print fx, fy
```

The condition is that at the call of the procedure the force array is initialized, either to zero or the gravitational forces only.

<hr>The addition of the term `eps**2`

in the computation of the distance serves to avoid zero distances and thus infinite or very large values in the force computation. This epsilon `eps`

should of course be small against typical distances of particles resp. the dimension of the scene, here `sigma=1.01`

is a typical distance, thus `eps=1e-6`

is small against it.