C++ : How do I make this particle simulation faster?

Question

I have a simple 2D simulation which has particles travelling at given velocities and could collide with each other and the walls. But the fps drops when there are more than 500 particles. I want to simulate at least 5000 particles smoothly because I'll be adding more features to the program. Is there a different and efficient approach to this? Thanks!

gas.cpp:

#include<cstdlib>
#include<vector>
#include<glut.h>
#include "gas.h"
#include<time.h>


float t = 1;        //time step

void gas::p(float pos_x, float pos_y, float vx, float vy, float mass)
{
    srand(time(0));
    m = mass;
    x = pos_x;
    y = pos_y;
    velx = vx;
    vely = vy;
    size = 3;
}


void gas::draw()
{
    glColor3f(1, 1, 1);
    glVertex2f(x, y);

}
void gas::move(float t)
{
    x += velx * t;
    y += vely * t;
}

float temp;

//Function to be ran at every frame:
void run(std::vector <gas>& bn)
{
    int it = 0;
    for (gas& i : bn)
    {
        int jt = 0;
        for (gas& j : bn)
        {

            if (it != jt)
            {   
                //Handling collisions:
                if (abs(i.y - (j.y + j.size)) < 1 && (abs(j.x - i.x) <= i.size + 1) && i.vely < 0)
                {
                    temp = i.vely;
                    i.vely = j.vely;
                    j.vely = temp;

                }
                if (abs(j.y - (i.y + i.size)) < 1 && (abs(i.x - j.x) <= j.size + 1) && i.vely > 0)
                {
                    temp = i.vely;
                    i.vely = j.vely;
                    j.vely = temp;  
                }
                if (abs(j.x - (i.x + i.size)) < 1 && (abs(i.y - j.y) <= i.size + 1) && i.velx > 0)
                {
                    temp = i.velx;
                    i.velx = j.velx;
                    j.velx = temp;
                    
                }
                if (abs(i.x - (j.x + j.size)) < 1 && (abs(j.y - i.y) <= i.size + 1) && i.velx < 0)
                {
                    temp = i.velx;
                    i.velx = j.velx;
                    j.velx = temp;
                }
            }

            jt += 1;
        }

        //Boundary Collisions:
        if (i.x > 600 - i.size) { i.x = 600 - i.size; i.velx = -i.velx; }
        if (i.x < i.size) { i.x = i.size; i.velx = -i.velx; }
        if (i.y > 600 - i.size) { i.y = 600 - i.size; i.vely = -i.vely; }
        if (i.y < i.size) { i.y = i.size; i.vely = -i.vely; }

        
        i.move(t);
        it += 1;

    }
}

gas.h:

#pragma once

class gas
{
public:

    float m = 1;
    float x = 0;
    float y = 0;
    float velx = 0;
    float vely = 0;
    float size = 3;
    float density = 100;
    float r = 1; float g = 1; float b = 1;


    void p(float pos_x, float pos_y, float vx, float vy, float mass);
    void draw();
    void move(float t);
};

void run(std::vector<gas>& bn);

simulate.cpp (main file):

#include<cstdlib>
#include<glut.h>
#include<vector>
#include "gas.h"
#include<thread>
#include<time.h>

void display();
void reshape(int, int);

const int n = 600;       //Number of particles

void init()
{
    glClearColor(0, 0, 0, 1);
}

std::vector <gas> b(n);

void timer(int)
{
    run(b);
    glutPostRedisplay();
    glutTimerFunc(1000 / 60, timer, 0);
}

void show(int argc, char** argv)
{
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowPosition(300, 60);
    glutInitWindowSize(600, 600);


    glutCreateWindow("Particles");
    glutDisplayFunc(display);
    glutReshapeFunc(reshape);

    glutTimerFunc(1000, timer, 0);      // Args: Time delay per frame in milliseconds, function to be called
            
    init();
    glutMainLoop();
}

int main(int argc, char** argv)
{   
    int it = 0;
    for (gas& i : b)
    {
        srand(it);
        i.p(rand() % 600, rand() % 600, (rand() % 10) * pow(-1, it + 1), (rand() % 10) * pow(-1, it), 1);
        it += 1;

    }
    show(argc, argv);
}


void display()
{
    glClear(GL_COLOR_BUFFER_BIT);
    glLoadIdentity();

    glPointSize(3);
    glBegin(GL_POINTS);

    for (gas& i : b)
    {
        i.draw();
    }
    glEnd();
    glutSwapBuffers();
}

void reshape(int w, int h)
{
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(0, 600, 0, 600);
    glMatrixMode(GL_MODELVIEW);
}