Starting velocities: Maxwell–Boltzmann

The idea

In MD, temperature is the average kinetic energy of the atoms. Hotter means faster. But at a given temperature the atoms do not all move at the same speed. Some are slow, some fast, and the speeds follow the Maxwell–Boltzmann distribution.

• Each velocity component (v_x, v_y, v_z) is a Gaussian centred on zero, with variance k_BT/m.
• The speed |v| then follows the skewed Maxwell–Boltzmann curve: no atoms at zero speed, a peak, and a long tail of fast atoms.
• Heavier atoms move slower at the same temperature: typical speed ∝ √(T/m).

The three-step recipe used by MD codes:

1. Draw each component from a Gaussian with variance k_BT/m.
2. Remove the drift: subtract the centre-of-mass velocity so the whole box does not fly off.
3. Rescale all velocities by √(T_target/T_now) so the temperature is exactly right.

Getting the exact shape right is less critical than it sounds. Collisions redistribute energy within a few picoseconds, so any reasonable start relaxes to Maxwell–Boltzmann. The Gaussian draw just gets you there immediately.

The algorithm

The recipe with masses, for N atoms:

for each atom i, each direction d:
    v[i, d] = sqrt(kB * T / m[i]) * gaussian()
v_cm = sum(m[i] * v[i]) / sum(m)        # mass-weighted
v -= v_cm                               # no net momentum
T_now = sum(m[i] * |v[i]|^2) / ((3N - 3) * kB)
v *= sqrt(T_target / T_now)

1. Use the mass-weighted centre-of-mass velocity when atoms have different masses.
2. Count 3N − 3 degrees of freedom after removing the centre-of-mass motion. Subtract more if bonds are constrained (Lesson 4.4).
3. Use a new random seed for each independent run. Different seeds give different trajectories, which is how you estimate error bars.

In LAMMPS this is velocity all create 300.0 4928459 dist gaussian mom yes. In GROMACS it is gen_vel = yes, gen_temp = 300.