Optimizing Monte Carlo simulation of particle movement in 2D space with C++

I am working on a program that simulates the movement of particles in a 2D space using a Monte Carlo method. The program should simulate the movement of each particle according to the following rules:

1. At each time step, each particle has a probability p of moving to a new location chosen randomly within a certain radius r of its current location.
2. If a particle moves to a location occupied by another particle, the two particles should "collide" and move to a new location chosen randomly within a certain radius s of their previous location.
3. If a particle reaches the edge of the 2D space, it should bounce off the wall and continue moving in the opposite direction.

I am having trouble implementing these rules in a way that is both efficient and accurate. Here is the code I have so far:

void move_particles(vector<Particle> &particles, double p, double r, double s) {
  for (int i = 0; i < particles.size(); i++) {
    if (should_move(p)) {
      // choose new location within radius r
      double new_x = particles[i].x + r * (2 * rand() - 1);
      double new_y = particles[i].y + r * (2 * rand() - 1);

      // handle boundary conditions
      if (new_x < 0) {
        new_x = -new_x;
      }
      if (new_x > WIDTH) {
        new_x = 2 * WIDTH - new_x;
      }
      if (new_y < 0) {
        new_y = -new_y;
      }
      if (new_y > HEIGHT) {
        new_y = 2 * HEIGHT - new_y;
      }

      // check for collisions with other particles
      for (int j = 0; j < particles.size(); j++) {
        if (i == j) continue;
        if (distance(new_x, new_y, particles[j].x, particles[j].y) < 2 * s) {
          // choose new location within radius s
          new_x = particles[j].x + s * (2 * rand() - 1);
          new_y = particles[j].y + s * (2 * rand() - 1);
        }
      }

      particles[i].x = new_x;
      particles[i].y = new_y;
    }
  }
}

This code seems to work for small numbers of particles, but when I try to run it with a large number of particles (e.g. 1000), the program becomes very slow and sometimes crashes. Can anyone help me understand why this might be happening and suggest ways to optimize the performance of the program?

>Solution :

The biggest performance issue is every loop is making n^2 checks, optimising the existing code will improve that but the biggest improvement would come from not checking every pair of particles. https://gameprogrammingpatterns.com/spatial-partition.html covers many options.

As for accuracy, you mention that particles should move to a new location on a collision, but you don’t check for subsequent collisions after moving a particle on its first collision, is that what you intended? (new_x and new_y may now collide with a previously checked particle)

As for your crash I dont see any problems in your code, where does your debugger say the crash is happening?