How can I optimize this compute shader?

I'm coding an N-body simulation in OpenGL.

This is the best performance I could achieve in my GTX960 4GB (frame rate limited at 60fps) and using the compute shader I provided here.

Num particles = 15.000; ms = 16.6ms; fps = 60fps
Num particles = 17.000; ms = 16.6ms; fps = 60fps
Num particles = 18.000; ms = 20.0ms; fps = 50fps
Num particles = 19.000; ms = 30.0ms; fps = 30fps
Num particles = 20.000; ms = 40.0ms; fps = 25fps
Num particles = 30.000; ms = 80.0ms; fps = 12.5fps

My question is if there's something I could do to optimize the simulation. Maybe by separating the struct and passing two arrays (one of pos and the other of velocities)?

I know that Barnes Hut is faster than this but I'm not sure if I can do recursion in compute shaders, if not then it would complicate a lot the implementation.

I was thinking to use a less complex optimization using uniform grids inspired by the grid collison detection algorithm. I don't know if it's a good idea.

#version 440
layout( local_size_x = 64, local_size_y =1, local_size_z = 1  ) in;


const float G = 6.673e-11;
const float SOFTENING = 9000.f;   // softening parameter (just to avoid infinities when calculating the force of gravity)

uniform int numParticles;
uniform float deltaTime = 1.f/60.f;

struct Particle{
    vec3 pos;
    vec3 velocity;
};


layout(std430, binding=0) buffer particlesBuffer {
    Particle particles[];
};


void updatePos(inout vec3 acceleration, inout int index){
    particles[index].velocity += acceleration * deltaTime * 0.1f;
    particles[index].pos += particles[index].velocity * deltaTime * 0.1f;
}


vec3 computeGravityForces(inout int index){
    vec3 total_force = vec3(0.0f, 0.0f, 0.f);

    // For each neighbor particle
    for(int j = 0; j < numParticles; j++){
        if(index != j){


            vec3 vec_i_j = particles[j].pos - particles[index].pos;
            vec3 normalized_vector = normalize(vec_i_j);

            const float squaredDistance = dot(vec_i_j, vec_i_j);
            const float strength = (G * (5.972f * 10e14)) / (squaredDistance + SOFTENING);

            // Add current force
            total_force += normalized_vector * strength;
        }
    }
    return total_force;
}


void updateParticle(inout int index){
    vec3 acceleration = computeGravityForces(index);
    updatePos(acceleration, index);
}

void main() {
    uint index = gl_GlobalInvocationID.x;
    updateParticle(index);

}