The code doesn't work properly because the simulation is blocked by an error:
ValueError: cannot convert float NaN to integer
So, I putted some print to view where the problem was and it happens in the x2, and y2 coordinates. From time to time they became NaN and the following assegnation crash: r = int(sqrt(((x2 - x1)**2) + ((y2 - y1)**2))) #Distance between the particles
I don't really understand where they are coming from.