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.

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. This would not, however, work as well if the computations are not so cleanly separable.

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.

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.