Velocity Maxwell-Boltzmann distribution for dummies I have a volume with $N$ molecules; I need to assign to each particle a velocity vector:
$$|\mathbf{v}_{i}|=[v_{x}, v_{y}, v_{z}]^{T}$$ for the $i$-th molecule; the velocities must follow the Maxwell-Boltzmann distribution; here is the question: I know the definition of the M.B. distribution, how do I extract the velocities that I have to assign to each molecule?
Basically what I want to know is how to calculate
$$v_{x}=?$$ $$v_{y}=?$$ $$v_{z}=?$$
for each molecule.
 A: The Maxwell-Boltzmann distribution is a continuous probability distribution of particles as a function of either momentum, energy, or speed. Since you want the velocity vector, we start with that:
$$
p_V(v_x,v_y,v_z)= \left(\frac{m}{2\pi kT}\right)^{3/2}\exp\left(-\frac{m}{2kT} v^2\right)
$$
where $v^2=v_x^2+v_y^2+v_z^2$. Note that, due to the product rule of exponential functions, you can actually write the above as
$$
p_V(v_x,v_y,v_z)=p_V(v_x)p_V(v_y)p_V(v_z)
$$
with
$$
p_V(v_x)=\left(\frac m{2\pi kT}\right)^{1/2}\exp\left[-\frac{m}{2kT}v_x^2\right]\tag{1}
$$
and likewise for $p_V(v_y)$ and $p_V(v_z)$. Thus, you can find the probability distribution of the vector velocity by finding three individual probabilities and multiplying them together.
In order to sample from the distribution in (1), you need the cumulative distribution function, defined as (for the continuous case)
$$
F_X(x)\dot=\int_{-\infty}^xp_V(v_x)\,dv
$$
The final result is
$$
F_X(x)=\frac12\left[1+{\rm erf}\left(\frac{x}{a\sqrt{2}}\right)\right]\tag{2}
$$
where $a=\sqrt{kT/m}$ and ${\rm erf}(\eta)$ is the error function. Note that $F_X(x)\in\left[0,1\right]$.
It is from this that you can draw your random velocities to match the distribution. Pick a random number in the range 0 to 1, set it as $F_X(x)$ and find the $x$ such that the two sides of (1) are equal. There are a number of ways to do this, you'll have to find the one that works best for you.
A: Here's a couple of points to add to Kyle Kanos's full answer, which will give you a couple more ways of computing and recalling facts related to what you want.
The Maxwell Boltzmann distribution happens to be a Chi distribution for three degrees of freedom, which means it is the distribution of $\sqrt{v_1^2+v_2^2+v_3^2}$, where the $v_i$ are all identical variance, zero mean Gaussian random variables. This in turn is related to the Chi-squared distribution, the probability density of the square of this quantity, i.e. $v_1^2+v_2^2+v_3^2$. 
So, to get the $v_i$, you need to make three independent zero-mean random numbers normally distributed with standard deviation $\sigma=\sqrt{\frac{k\,T}{2\,m}}$. Two other methods, further to Kyle's manipulation of the cumulative distribution, are the following:
First Method, using the Box Muller transformation: The Box-Muller transformation as described by:
Everett Carter, "Generating Gaussian Random Numbers"
is an excellent implementation of this principle. 
Press, W.H., B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, 1986; Numerical Recipes, The Art of Scientific Computing, Cambridge University Press, Cambridge
also describes this in full. In short, we take two independent random numbers $x_1$ and $x_2$, which are uniformly distributed in the interval $[0,\,1]$. So you simply need a uniform random number generator. Then we form:
$$y_1 = \sqrt{\frac{k\,T}{2\,m}}\,\sqrt{-2\,\log{x_1}}\,\cos(2\,\pi\,x_2)$$
$$y_2 = \sqrt{\frac{k\,T}{2\,m}}\,\sqrt{-2\,\log{x_1}}\,\sin(2\,\pi\,x_2)$$
and then $y_1$ and $y_2$ are Gaussian variables with the required variance $\frac{k\,T}{2\,m}$ and mean nought. So you would need to call this method three times to get two velocity vectors (it would yield six independent Gaussian variables which you can assemble into the two vectors).
Second method: Use the central limit theorem. Add, say, one hundred independent random numbers, each uniformly distributed in $[0,\,1]$. The sum $\Sigma$ will be very nearly Gaussian with a mean of $50$ and a variance of $25/3$ (equal to 100/12, since the variance of the uniform variable is $1/12$). So now you compute:
$$\sqrt{\frac{3\,k\,T}{50\,m}} (\Sigma-50)$$
to get each velocity component.
A: I know I'm answering this question almost three years late but I've encountered the same problem and want to share my solution, in case anybody else finds this post, as I did. The previous answers give excellent background on the Maxwell-Boltzmann distribution (MBD) so I'm not going to go into detail on that.
Once you find your preferred method to obtain "random" speeds distributed according to the MBD (I got a clue from this video youtu.be/C82JyCmtKWg?t=9m4s) what you can do is generate two other random numbers. So, for each speed $v_i$, you also compute $\theta_i$ and $\phi_i$, that way you can do the following vector decomposition:
$$
v_{i,x} = v_i \sin(\phi)\cos(\theta) \\
v_{i,y} = v_i \sin(\phi)\sin(\theta)\\
v_{i,z} = v_i \cos(\phi)
$$
