Use `variable_timestep` in equilibration phase

**Is your feature request related to a problem? Please describe.**
One issue that comes with using the initial random parameters is that sometimes the parameters will cause the subject system to "explode" (i.e. move so fast that the integration algorithm destabilises and fundamental MD assumptions don't hold).

**Describe the solution you'd like**
One potential solution to this would be to use a variable time-step algorithm in the equilibration phase. This would allow the system to hopefully reach a stable state.

This could be implemented as part of @CosmicPidgeon 's recovery from failures PR (#1354 or more likely a subsequent one building on it). On failure, reducing the time-step may be a viable solution (with issued warning).

DLPoly also provides (and LAMMPS may too) a `variable_timestep` feature which computes relative to the maximum velocities/forces in the system a suitable timestep dynamically as part of the calculation. However, due to the necessary statistics, this should only be enabled during the equilibration phase or the statistics of the Observables will difficult to compute.