<p>I assume that you want to simulate a gas of particles with potential, which has attractive and hard-core repulsive parts.</p>

<p>Molecular dynamics simulation of many-particle systems is a well-developed area of science. The hard part of such simulation is to sample particles initially in such way that they don't overlap. This is solved in the classical paper by <a href="http://bayes.wustl.edu/Manual/EquationOfState.pdf" rel="nofollow noreferrer">Metropolis et al</a>. But you neglect this part anyway. There are also many textbooks and papers explaining, how to perform molecular dynamics simulations with hard spheres, like <a href="https://arxiv.org/pdf/1211.6718.pdf" rel="nofollow noreferrer">this</a> or <a href="https://www.google.com/url?sa=t&amp;rct=j&amp;q=&amp;esrc=s&amp;source=web&amp;cd=2&amp;ved=0ahUKEwis99Oo7OnXAhVVzGMKHbC6DCMQFggvMAE&amp;url=https%3A%2F%2Fdspace.library.uu.nl%2Fbitstream%2Fhandle%2F1874%2F8984%2Ffrenkel_89_molecular_hard_particles.pdf%3Fsequence%3D1&amp;usg=AOvVaw30eGrhRL6OgtoDLHSb1Qn_" rel="nofollow noreferrer">this</a> or the <a href="http://aip.scitation.org/doi/pdf/10.1063/1.1730376" rel="nofollow noreferrer">1959 paper, where they first used Molecular Dynamics</a>.</p>

<p>If your project is to show nice spheres that collide like in Windows screensaver, just make your hard-core potential softer, decrease your time step and it will work like a charm. If you aim at scientific rigour, read those papers and don't forget to control energy conservation. Soft repulsion can be added, for example, by replacing</p>

<pre><code>double S_tr = 1 /r_ij_square;
</code></pre>

<p>by</p>

<pre><code>double S_tr = 1.0 / r_ij_square - 1.0 / (r_ij_square * r_ij_square);
</code></pre>

<p>The code formatting, structure and variable naming leave much to be desired. I would recommend to try <a href="https://clang.llvm.org/docs/ClangFormat.html" rel="nofollow noreferrer">clang-format</a> tool. Structurewise, there is a lot of duplication of periodic boundary conditions.</p>

<p>With the repulsive potential you may get a system that experiences a phase transition, an exciting object to study. Good luck!   </p>

<p>UPDATE:</p>

<p>Note that the following piece of code is bugged:</p>

<pre><code>double rho, v0, tr, xr, l0, eta, dt, x[N], y[N];
double L = pow(N / rho , 0.5), l0_two = l0 * l0;
// ...
l0 = 1;
</code></pre>

<p>You first compute l0_two = l0 * l0 (never start variable name with l, it's really easy to confuse it with 1, sqrt(x) should be preferred to pow(x, 0.5)) and then assign l0 = 1. Therefore l0_two is undefined and is probably assigned to 0. This may cause part of your problem.</p>

<p>After fixing that bug you can go for two solutions. The first would be the soft repulsive potential described above if particles are close enough. The second is to do exactly as in the <a href="https://arxiv.org/pdf/1211.6718.pdf" rel="nofollow noreferrer">paper with event-based MD cited above</a>.</p>
How to add hard-sphere repulsion in my simulation of attractive particles?
<p>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.</p>

<p>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?</p>

<p>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.</p>

<p>This is my code:</p>

<pre><code>#include &lt;iostream&gt;
#include &lt;math.h&gt;
#include &lt;vector&gt;
#include &lt;array&gt;
#include &lt;list&gt;
#include &lt;random&gt;
#include &lt;functional&gt;
#include &lt;fstream&gt;
#include &lt;string&gt;
#include &lt;sstream&gt;
#include &lt;algorithm&gt;
#include &lt;chrono&gt;
#include &lt;set&gt;

using namespace std;

    std::minstd_rand gen(std::random_device{}());
    std::uniform_real_distribution&lt;double&gt; unirnd(0, 1);

double PBCtwo(double pos, double L)
    {
    if(pos &gt; L / 2.0)
        return pos-L;
    else if (pos &lt; -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 &lt;&lt; to_string(N) &lt;&lt; "\n";
    configFile &lt;&lt; to_string(L) &lt;&lt; "\n";
    for (int i = 0; i &lt; N; i++)
    {   x[i] = unirnd(gen) * L;
        y[i] = unirnd(gen) * L;
        configFile &lt;&lt; to_string(x[i]) &lt;&lt; "\t" &lt;&lt; to_string(y[i]) &lt;&lt;   "\n";
    }
    configFile.close();

    while (c &lt; c_prod)
        {
        double dx[N], dy[N];
        c++;
        for(int i = 0; i &lt; N; i++)
           {
            dx[i] = 0;
            dy[i] = 0;
            double S_try = 0.0, S_trx = 0.0;
            for(int j = 0; j &lt; 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 &gt; 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 &lt; N; i++)
            {
            x[i]+=  dt * dx[i];
            y[i]+=  dt * dy[i];
            if (x[i] &gt; L){
              x[i]-= L;} 
            else if( x[i] &lt; 0) {
            x[i]+= L;}
            if (y[i] &gt; L){
             y[i]-= L;} 
            else if( y[i] &lt; 0){
                y[i]+= L;}
            }
        }
    ofstream finalConfigFile;
    finalConfigFile.open ("Final Configuration.txt");
    finalConfigFile &lt;&lt; to_string(N) &lt;&lt; "\n";
    finalConfigFile &lt;&lt; to_string(L) &lt;&lt; "\n";
    for (int i = 0; i &lt; N; i++)
        {
        finalConfigFile &lt;&lt; to_string(x[i]) &lt;&lt; "\t" &lt;&lt; to_string(y[i]) &lt;&lt;"\n";   
        }
    finalConfigFile.close();
    return 0;
}
</code></pre>