如何在我对有吸引力的颗粒的模拟中添加硬球排斥

我模拟了一个互相吸引的粒子系统。吸引力的强度取决于颗粒之间的距离。边界条件和相互作用是周期性的。由于吸引力，粒子彼此相处并聚集了一个圆圈。

我想添加硬球的排斥，以便每当两个或多个粒子聚集在相同位置时，我希望它们在连接其中心的线路中分开，直到它们不重叠为止。我该怎么做？

当存在吸引人的相互作用时，添加硬球的情况比通常的情况要困难，因为在某些情况下可能有4个或更多粒子在同一位置。

这是我的代码：

#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;
}