Anne M. Brant I think sampling the velocities both angular and linear using Maxwell-Boltzmann distribution at the desired temperature might be a solution for your problem.

If all the bonds in your molecule have fixed lengths and angles, then you can write down a fixed inertia tensor for your system. Identify axes of rotation about the center of mass that diagonalize your inertia tensor, and then draw from a Gaussian distribution for the linear or angular velocity of each translational and rotational coordinate. Scale this based on the temperature and the mass or moment corresponding to a given coordinate. Regardless of species, you should have an expected kinetic energy of 1/2 k_BT associated with each degree of freedom.

If you instead have a more complicated system- perhaps your bond lengths are fixed, but some bond angles, or maybe dihedral angles, can vary- I've found this paper helpful for generating a thermodynamically accurate distribution of translational and rotational velocities for arbitrary molecules:

Jain, Abhinandan, In-Hee Park, and Nagarajan Vaidehi. "Equipartition principle for internal coordinate molecular dynamics." Journal of chemical theory and computation 8.8 (2012): 2581-2587.

I would make sure that you have 1/2kT of kinetic energy per degree of freedom on each molecule. If it is a rigid body simulation, you have 3 translational DoFs and 3 rotational DoFs, find the required velocity and angular velocity for each molecule and generate randomly for each molecule. If you want your center of mass not to be moving, you have to remove the total linear velocity of the system from each molecule and then rescale the velocities to get your total required kinetic energy back. Also, rescale all the velocities to get you traslational kinetic energy by (N - 3)/N to reflect that you have 3 DoFs less due to the resting center of mass.

I would not bother to sample the velocities from a distribution, though, they should get there on their own after a short simulation.