In this article, we investigate molecular dynamics (MD) trajectories of a butane molecule, as obtained using different types of thermostats. Results show that at low temperature, where the harmonic approximation holds, the Nose'−Hoover (NH) thermostat fails to reproduce the statistical mechanical behavior, even using simulation lengths of millions of time steps, whereas the Gaussian isokinetic (IG) thermostat reproduces quite well the expected statistical mechanical values. The Berendsen's coupling (BC) provides good results for basic properties such as the average potential and kinetic energies but fails in reproducing the canonical fluctuations. Moreover, using the speed of divergence of initially nearby trajectories in phase space as a measure of the dynamical chaoticity, we found that the NH thermostat provides very slow divergence for the physical phase space degrees of freedom, concentrating most of its chaoticity in the dynamics of the thermostat virtual degree of freedom. On the contrary, the IG thermostat provides always highly diverging trajectories in phase space, characterized by a high chaoticity of each degree of freedom. Finally, the BC thermostat provides a moderate chaotic behavior for all of the degrees of freedom. Such results suggest that even assuming for both the “rigorous” algorithms (NH and IG) a full ergodic behavior the NH thermostat could require an extremely long time to achieve convergence of the time averaged properties.