Active Brownian Motion I am attempting to write a Python code to simulate many particles in a confined box. These particles behave in such a way that they move in the box in straight lines with a slight angular noise (small changes in the direction of the particle path). They should interact by acknowledging the other particle and 'shuffle/squeeze' past each other and continue on their intended path, much like humans on a busy street. Eventually, the particles should cluster together when the density of particles (or packing fraction) reaches a certain value.
However, I have a feeling there are parts of the code that are inefficient or which could be either sped up or written more conveniently.
If anyone has any improvements for the code speed or ideas which may help with the interactions and/or angular noise that would be much appreciated. I will also leave an example of an animation which is my aim: https://warwick.ac.uk/fac/sci/physics/staff/research/cwhitfield/abpsimulations
The above link shows the animation I am looking for, although I don't need the sliders, just the box, and moving particles. The whole code is shown below:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


def set_initial_coordinates():
    x_co = [np.random.uniform(0, 2) for i in range(n_particles)]
    y_co = [np.random.uniform(0, 2) for i in range(n_particles)]
    return x_co, y_co


def set_initial_velocities():
    x_vel = np.array([np.random.uniform(-1, 1) for i in range(n_particles)])
    y_vel = np.array([np.random.uniform(-1, 1) for i in range(n_particles)])
    return x_vel, y_vel


def init():
    ax.set_xlim(-0.05, 2.05)
    ax.set_ylim(-0.07, 2.07)
    return ln,


def update(dt):
    xdata = initialx + vx * dt
    ydata = initialy + vy * dt
    fx = np.abs((xdata + 2) % 4 - 2)
    fy = np.abs((ydata + 2) % 4 - 2)

    for i in range(n_particles):
        for j in range(n_particles):
            if i == j:
                continue

            dx = fx[j] - fx[i]  # distance in x direction
            dy = fy[j] - fy[i]  # distance in y direction
            dr = np.sqrt((dx ** 2) + (dy ** 2))  # distance between x

            if dr &lt;= r:
                force = k * ((2 * r) - dr)  # size of the force if distance is less than or equal to radius

                # Imagine a unit vector going from i to j
                x_comp = dx / dr  # x component of force
                y_comp = dy / dr  # y component of force

                fx[i] += -x_comp * force  # x force
                fy[i] += -y_comp * force  # y force
    ln.set_data(fx, fy)
    return ln,


# theta = np.random.uniform(0, 2) for i in range(n_particles)
n_particles = 10
initialx, initialy = set_initial_coordinates()
vx, vy = set_initial_velocities()
fig, ax = plt.subplots()
x_co, y_co = [], []
ln, = plt.plot([], [], 'bo', markersize=15)  # radius 0.05
plt.xlim(0, 2)
plt.ylim(0, 2)

k = 1
r = 0.1
t = np.linspace(0, 10, 1000)

ani = FuncAnimation(fig, update, t, init_func=init, blit=True, repeat=False)
plt.show()
Coordinates as a matrix
This:
x_co = [np.random.uniform(0, 2) for i in range(n_particles)]
y_co = [np.random.uniform(0, 2) for i in range(n_particles)]
return x_co, y_co

is a little baffling. You call into Numpy to generate uniformly distributed numbers (good) but then put it in a list comprehension (not good) and separate x and y into separate lists (not good).
Consider instead
return np.random.uniform(low=0, high=2, size=(2, n_particles))

and similarly for set_initial_velocities:
return np.random.uniform(low=-1, high=1, size=(2, n_particles))

All of the rest of your code, such as the calculations for xdata/ydata should similarly vectorize over the new length-2 dimension eliminating distinct code for x and y.