Unfortunately, for a many-particle system like BzAr 
rather long
trajectories are needed to estimate thermodynamic properties of the
microcanonical ensemble (constant E, V and N) using
Molecular Dynamics (MD) with any accuracy. As for Monte Carlo studies
simulations are likely to be trapped in isolated regions of the constant energy
hyperplane in phase space.
However, insurmountable transition states can make simulations restricted to
such regions meaningful, but even calculating locally ergodic trajectories can
be time consuming.
Our MD integration in ORIENT 3.2 [27] is done with an adaptive step size propagator due to Burlish and Stoer [38, 39]. Here the result for a large integration step is extrapolated from a sequence of finer and finer substep-integrations of the large interval. The extrapolation procedure provides us then with an error estimator for each integration step. This is used for the dynamical adjustment of the step size. Both the chosen propagator and the adaptive step size are designed to minimise the number of force calculations. This strategy is much more efficient then corrector-predictor and static step size schemes. In addition it allows us to directly control the integration quality with regard to the integration variables and the integrals of motion. In our adaptive step size control we consider both the accuracy of the phase space coordinates and the conservation of energy.
In the microcanonical simulations we loop over increasing values of the total energy. At each new energy we start from the last configuration of the previous energy. The initial kinetic energies are equally distributed amongst the vibrational degrees of freedom only but spread out into the internal rotations after a few time-steps. The initial centre of mass motion and the overall angular momentum are always corrected so that they vanish at each new energy. These overall momenta are well conserved in the integration since we do not use a constraining container for MD.
All trajectories are monitored and quenched periodically to provide insight
into the region of phase space in question. We calculate
the kinetic temperature 
as,
resulting from the average kinetic energy per degree of freedom [40]. Here conservation of the overall linear and rotational momenta has to be taken into account.
We treat the benzene molecule as a rigid body in all simulations. The orientation space is represented by a Euclidean metric using the corresponding quaternion parameters [37, 41, 42, 43]. An equivalent angle-axis description [43] is also possible without singularities, in contrast to e.g. Euler angles. The latter approach is convenient for sampling orientation space in Monte Carlo simulations and we used it in static step size calculations. However, the integration scheme with adaptive step size control is much more convenient in the quaternion representation where larger steps can be taken without problems due to non-commuting finite rotations.
