Simulation of interacting particles I'm simulating a system of particles which interact with each other and move and rotate due to the interaction.  There are 2000 particles in total but each particle interacts with other particles and their image (I'm using a periodic boundary condition) if particles and their images have distances larger thanl0. I just consider image of particles in 8 near boxes to the main box.

y[N],x[N] are the position of points in the plane; dx[n] and dy[N] their changes in each time step. c shows different time steps. e_x and e_y shows x and y component of direction of each particle respect to x and y-axis.

int main()
{
  unsigned short int N, l0 , R0;
  double rho, v0, eta, dt, a, b, D;
  long c_equilib, c_prod, c_display;
     while (c &lt; 1000000)
            {
         c++;
        for(unsigned short int i=0; i&lt;2000; i++){
          double S_theta=0.0, S_x=0.0, S_y=0.0, S;
                    for(unsigned short int j=0; j&lt;2000; j++)
                    {
                        if (j==i) continue;
                        for(int ky=-1; ky&lt;=1; ky++)
                        {
                            for(int kx=-1; kx&lt;=1; kx++)
                            {
                                double r_x_ij=x[i]-(x[j]+kx*L),
                                    r_y_ij=y[i]-(y[j]+ky*L),
                                r_ij=sqrt(pow(r_x_ij,2.0)+pow(r_y_ij,2.0));
                                if (r_ij &gt;= l0)
                                { 
                                    r_x_ij/=r_ij;
                                    r_y_ij/=r_ij;
                                    S_theta += (e_x[i]*e_x[j]+e_y[i]*e_y[j]) / pow(r_ij, 3) * ( 3.0* (r_x_ij*e_x[j]+r_y_ij*e_y[j]) *
                                                            pow(e_x[i]*r_y_ij - e_y[i]*r_x_ij, 2.0) -
                                                            e_x[i]*e_y[j] + e_y[i]*e_x[j] );

                                    S = ( 3.0* pow(r_x_ij*e_x[j]+r_y_ij*e_y[j], 2.0) -1.0 ) / pow(r_ij, 2.0);
                                    S_x += r_x_ij*S;
                                    S_y += r_y_ij*S;
                                }
                            }
                        }
                    }
        }
dx[i] += dt*a*S_x;
            dy[i] += dt*a*S_y;
        }
e_x[i]=cos(theta[i]);
e_y[i]=sin(theta[i]);
x[i]+=dx[i];
y[i]+=dy[i];
}
}


I'm looking to optimize the large loop of this C++ program.  This is one part of my simulation (I have ignored some calculations inside the loop to make the question simpler). This part is extremely time-consuming. How can I make it faster?
I don't know how fast you expect this to run, but you need to understand that the inner part of your nested loop runs 36 trillion (1000000 * 2000 * 2000 * 3 * 3) times. This means that small changes in the run time of the inner loop will have a large impact on the total run time of the code, turning 1 hour of runtime to years.

Without knowing anything about the domain of your simulation, I can only propose this:


Reduce the number of times the loops run
Parallelize the code, so that you can have multiple threads working on  different parts concurrently.
Optimize the inner loop as much as possible. Perhaps there is some approximation of the calculation you're running that's significantly faster? Is there any way that you can make better use of processor caches or avoid branch mispredictions? Have a look at this latency chart.