How to add hard-sphere repulsion in my simulation of attractive particles?

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 particles gather in the same position, I want them to seperate in the line linking their centers, till they don't overlap. How can I do this?

The situation for adding hard-sphere when there is an attracting interactions is harder than the usual case, as there could be some situations in which 4 or more particles 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;
}