Nested loop over std::vector

In my simulation I have a particle and a nbody class.

struct Particle{
    double m;           // mass
    double x[DIM];      // position
    double v[DIM];      // velocity
    double F[DIM];      // force 
    double F_old[DIM];  // force past time step
    void update_position(double dt);
    void update_velocity(double dt);
    void update_force(Particle*, Particle*);
};

class Nbody {
    private:
        const double dt;                // step size
        const double t_max;             // max simulation time
        std::vector<Particle> p;
    public:
        Nbody(unsigned int n_, double dt_, double t_max_);
        void comp_force();
        void comp_position();
        void comp_velocity();
};

I update the position of every particle by looping over every particle and calling void Particle::update_position().

void Nbody::comp_position() {
    for (Particle& particle: p) { 
        particle.update_position(dt);
    }
}

In the following code snippet, I show the update rule to compute the new position.

void Particle::update_position(double dt) {
    const double a = dt * 0.5 / m;
    for (unsigned int d=0; d<DIM; ++d) {
        x[d] += dt * (v[d] + a * F[d]);
        F_old[d] = F[d];
    }
}

However, I have problems with the force calculation where I want to compute the interaction of particles and therefore need a nested loop. How can I do that? Here is my approach so far, which gives me an error.

void Particle::update_force(Particle* i, Particle* j) {
    double r=EPS; // smoothing
    for (unsigned int d=0; d<DIM; d++) {
        r += sqr(j->x[d] - i->x[d]);
    }
    const double f = i->m * j->m / (sqrt(r) * r);
    for (unsigned int d=0; d<DIM; d++) {
        i->F[d] += f * (j->x[d] - i->x[d]);
        j->F[d] -= f * (j->x[d] - i->x[d]);
    }
}

void Nbody::comp_force() {
    for (Particle& particle: p) {
        for (unsigned int d=0; d<DIM; ++d) {
            particle.F[d] = 0.;
        }
    }
    for (unsigned int i=0; i<n; i++) {
        for (unsigned int j=i+1; j<n; j++) {
            update_force(&p[i], &p[j]);
        }
    }
}

Can I write this similar as the update_position method?