I have simulated a 2d system of particles which attract each other. Strength of attraction depends on distance between particles. The boundary condition and interactions are periodic. Because of attraction, particles go to each other and gather in a circle. I want to add hard-sphere repulsion so that, whenever two or more circles(particles gather in the same position, I want them to seperate in the line linking their centers till they don’t over lap.). consider raduce of each particle is constant called ‘a’. How can I add hard sphere repulsion between spheres? (Whenever PARTICLES overlap during the simulations, I separate them along the line connect- ing their centers.)The situation for adding hardsphere when there is an attracting interactions is harder than the usual case, as there could be some situations in which 4 or more particles come in the same position. This is my code:

#include <iostream> #include <math.h> #include <vector> #include <array> #include <list> #include <random> #include <functional> #include <fstream> #include <string> #include <sstream> #include <algorithm> #include <chrono> #include <set> using namespace std;     std::minstd_rand gen(std::random_device{}());     std::uniform_real_distribution<double> unirnd(0, 1); double PBCtwo(double pos, double L)     {     if(pos > L / 2.0)         return pos-L;     else if (pos < -L /2.0)         return L + pos;     else         return pos;     }  // main function int main() {   long c = 0;     int N=4000;      double rho, v0, tr,xr,l0, eta,dt, x[N],y[N],L=pow(N / rho , 0.5),l0_two = l0 * l0;                 rho = 2;                 v0 =300;eta = 1;dt = 0.0001;l0 = 1; c_prod = 500;c_display = 100;tr = -0.4;             // write initial configuration to the file     ofstream configFile;     configFile.open ("Initial Configuration.txt");     configFile << to_string(N) << "\n";     configFile << to_string(L) << "\n";     for (int i = 0; i < N; i++)     {   x[i] = unirnd(gen) * L;         y[i] = unirnd(gen) * L;         configFile << to_string(x[i]) << "\t" << to_string(y[i]) <<   "\n";     }     configFile.close();      while (c < c_prod)         {         double dx[N], dy[N];         c++;         for(int i = 0; i < N; i++)            {             dx[i] = 0;             dy[i] = 0;             double S_try = 0.0, S_trx = 0.0;             for(int j = 0; j < N; j++)               {                 if (j==i) continue;                    double delta_x = x[i]-x[j],                             delta_y = y[i]-y[j];                        double r_x_ij = PBCtwo(delta_x,L),                               r_y_ij = PBCtwo(delta_y,L),                               r_ij_square = r_x_ij * r_x_ij + r_y_ij * r_y_ij;                        if (r_ij_square > l0_two)                           {                            double r_ij = sqrt(r_ij_square);                           r_x_ij/= r_ij;                            r_y_ij/= r_ij;                           double S_tr = 1 /r_ij_square;                           S_trx += r_x_ij * S_tr;                           S_try += r_y_ij * S_tr;                           }                 }             dx[i] += tr * S_trx;               dy[i] +=  tr * S_try;             }         for(int i = 0; i < N; i++)             {             x[i]+=  dt * dx[i];             y[i]+=  dt * dy[i];             if (x[i] > L){               x[i]-= L;}              else if( x[i] < 0) {             x[i]+= L;}             if (y[i] > L){              y[i]-= L;}              else if( y[i] < 0){                 y[i]+= L;}             }         }     ofstream finalConfigFile;     finalConfigFile.open ("Final Configuration.txt");     finalConfigFile << to_string(N) << "\n";     finalConfigFile << to_string(L) << "\n";     for (int i = 0; i < N; i++)         {         finalConfigFile << to_string(x[i]) << "\t" << to_string(y[i]) <<"\n";            }     finalConfigFile.close();     return 0; }