What temperature means in MD

The idea

In a real gas, temperature measures how vigorously the molecules move. MD uses the same idea: the kinetic temperature is defined by the average kinetic energy. Every independent way an atom can move (a degree of freedom) carries on average ½k_BT of kinetic energy. This is equipartition. So we add up the kinetic energy and divide by the number of degrees of freedom.

Counting degrees of freedom. N atoms in 3D have 3N velocity components. Our simulations keep the total momentum at zero (Lesson 2.3), so 3 of them are not free: N_f = 3N − 3. Constraints remove more: a rigid water molecule has 6, not 9. Getting N_f wrong biases T, by 1% for 100 atoms and much more for constrained molecules.

T fluctuates. The kinetic energy at any instant is a sum over a finite number of atoms, so the instantaneous T jitters around its average. The relative size of the jitter shrinks as 1/√N. A 100-atom system flickers by several per cent; a million-atom system hardly at all. Only the time average is "the temperature".

T drifts at the start. If you start atoms on a perfect lattice (zero potential energy stored in displacements) with velocities at T₀, energy flows from kinetic into potential as the atoms leave their sites. In NVE the total energy is fixed, so T falls, often to about half of T₀ for a solid near equilibrium. This is why we equilibrate, and why we need thermostats (Lessons 3.3–3.5) to hold T at a target.

The algorithm

Measuring and setting the temperature:

def temperature(v, m):
    KE = 0.5 * sum(m * v**2)           # all atoms, all directions
    Nf = 3 * N - 3 - n_constraints
    return 2 * KE / (Nf * kB)

# set an exact starting temperature
v = gaussian(0, sqrt(kB * T0 / m))     # Maxwell-Boltzmann (Lesson 2.3)
v -= mean(v)                           # zero total momentum
v *= sqrt(T0 / temperature(v))         # rescale to exactly T0

1. Use the same N_f everywhere. MD codes print it (GROMACS reports the "number of degrees of freedom"; LAMMPS: compute_modify thermo_temp dof).
2. Rescaling all velocities by λ multiplies T by λ². This is the simplest thermostat, and the starting point for Berendsen (Lesson 3.3). In NVE, part of the added energy then flows into potential energy, so T settles below λ²T (in a solid, about halfway).
3. Do not judge T from one frame. Average over many steps, after the initial drift has gone.
4. After a lattice start, either start at about 2T_target or, better, equilibrate with a thermostat before the production run.