double sum_sin = 0.0, sum_cos = 0.0;
            int count = 0;
    #pragma omp parallel for reduction(+ : count ,sum_sin,sum_cos)
            for (vector<int>::iterator it = box_neighbors[bx[i]].begin(); it != box_neighbors[bx[i]].end(); ++it)
                {
                for (vector<int>::iterator itp = box_particles[*it].begin(); itp != box_particles[*it].end(); ++itp)
                     {
                     if(dist(x[i], y[i], x[*itp], y[*itp], L) < R0_two)
                        {
                        sum_sin+= sin(theta[*itp]);
                        sum_cos+= cos(theta[*itp]);
                        count+=1;  //number of neighbours of i'th particle
                        }
                     }
                }
            sum_sin/= count;
            sum_cos/= count;