In addition to G. Sliepen’s idea, you could use the STL’s `std::valarray<double>`.  This would let you replace something like

    for (int d = 0; d < DIM; ++d) {
        p->x[d] += dt * (p->v[d] + a * p->F[d]);
        p->F_old[d] = p->F[d];
    }

with something like

    p->F_old = p->F;
    p->x += dt * (p->v + a * p->F);

It would also be possible to lay out a structure of arrays rather than an array of structures.  If there are more particles than dimensions, this could let you perform wider vector operations on all the x-coordinates, then all the y-coordinates and all the z-coordinates, rather than being limited to the width of the coordinate system.  That is, each `p` might have only two or three parallel computations, but if you have a number of `std::array<std::valarray<double>, DIM>` with the x-coordinates in `x[0]`, the y-coordinates in `x[1]` and the z-coordinates in `x[2]`, the velocities in `v[0]`, etc., that might look like:

    for (size_t i = 0; i < x.size(); ++i) {
      F_old[i] = F[i];
      x[i] += dt * (v[i] + a * F[i]);
    }

and be able to use the full width of your vector registers.