I don't quite understand this suggestion. VelocityBoltzmann can already deal with multidimensional systems (aribtrary ndofs and natoms). Are you suggesting that it automatically restricts the velocity distributions that are generated to a subset of atoms?

I am asking as I am currently working on that code.

My idea was to make generating a Boltzmann distribution more convenient when a NQCDynamics.Simulation is already initialised.
So instead of VelocityBoltzmann(T, sim.atoms.masses, size(sim)), enable VelocityBoltzmann(T, sim) to generate the same initial distribution.

This makes sense. The solution I came up with is also not particularly nice. it works for now, but the interface could be easily improved.

This needs to be implemented in NQCDynamics.jl unless we want to draw Simulations out into a separate package.

Would also be useful to account for frozen atoms in simulations automatically → Generate guaranteed 0 velocity distributions for those atoms instead of Boltzmann. #12 makes multiple types of distributions for a UnivariateArray possible, so this can be done automatically in NQCDynamics.jl or manually by starting with a Matrix{UnivariateDistribution}.