Randomly initializing ensemble of particles in "toy" computer model I'm programming a "toy" model, and want to intialize the $(\mbox{positions},\mbox{momenta})$ of an ensemble of particles in three-space, using a uniform (pseudo) random number generator. But I'm a little lost about that: suppose the ensemble's characterized by a macroscopic thermodynamic temperature $T$. Then how do you use a uniform rng to realistically (reasonably realistically) distribute particle kinetic energies that reproduce the corresponding partition function?
Also, suppose the ensemble's constrained inside a finite box. I'd intended to just uniformly distribute initial positions randomly in $(x,y,z)$. But do the box boundaries affect that, or not? Should the (average) number-density of particles be a function of distance from the walls, or what? And ditto for momentum direction, i.e., I'd intended to use https://math.stackexchange.com/questions/44689/ to uniformly distribute initial directions in each solid angle. How (un)realistic is that, and what's the right way to account for the box boundary?
 A: How to initialise the velocities?
The way that molecular dynamics codes do this is to randomly draw numbers from the Maxwell-Boltzmann distribution, i.e. each of the three velocity components are drawn from the distribution
$$ p(v)=(2\pi mk_BT)^{-3/2}\exp\left(-\frac{mv^2}{2k_BT}\right) $$
Note that you can generate Gaussian-distributed random numbers using the Box-Muller transform.
How to initialise the positions?
If the box boundaries are periodic or reflective then they have no effect on the probability distribution, so a uniform distribution would be fine.
Although if your particles interact then you must be careful to avoid overlaps (and you should equilibrate before any ensemble sampling).
A: The general rule is that the distribution in phase space is weighted by $e^{-E/kT}$ where $E$ is the energy. 
In particular, if you're thinking about an ideal gas, there is no interaction between gas molecules and the wall (except at the moment the molecules bounce off), so the energy doesn't depend on position and the position distribution in the box is uniform. 
In a more realistic model of the wall, say one where the wall repels particles starting a distance $\lambda$ away, you would have a nontrivial potential $U(\mathbf{r})$, and draw positions according to $e^{-U(\mathbf{r})}$. But if you're using an ideal wall in your simulation the uniform distribution is perfectly correct.
Since kinetic energy is proportional to $p^2$, the momentum distribution goes as $e^{-p^2}$. Integrating over angles, you can choose the angle uniformly and randomly, then choose the magnitude, which is distributed as $p^2 e^{-p^2}$. 
However, choosing a random angle is a bit messy; an easier way is to note that $e^{-p^2}$ factorizes as
$$e^{-p^2} = e^{-p_x^2} e^{-p_y^2} e^{-p_z^2}$$
so you can just choose each component of the momentum independently, as a Gaussian. Again, since the kinetic energy doesn't depend on the position, you don't have to do any adjustment for 'being near the wall'.
