LAMMPS from scratch · 20 min

Velocities

One line, velocity all create 3.0 87287, gives 4000 atoms their starting velocities in four steps: random numbers, no net drift, optional no rotation, and an exact rescale to the target temperature. You read those steps in the source, check the kinetic energy by hand, and watch the velocity distribution change shape in the first 100 steps.

The four steps of velocity create, in velocity.cpp
Why the kinetic energy at step 0 is exactly (3N − 3)/2 × T, and what mom no does
Watch uniform random velocities become a Maxwell distribution through collisions

1 · Understand: four steps in one line

Step	What happens	Keyword (default)	velocity.cpp
1. random	each atom gets three random numbers, from a flat range −½…½ or a bell curve, divided by √m	`dist uniform` or `gaussian` (uniform)	lines 328–353
2. no drift	subtract the average velocity, so the whole crystal does not fly away	`mom yes` (yes)	line 357
3. no spin	remove any overall rotation	`rot yes` (no)	line 358
4. rescale	measure the temperature, multiply every velocity by √(T_target/T_measured)		lines 364–368

With loop geom (as in melt) the random numbers of each atom are seeded from the seed and the atom's own position (random->reset(seed,x[i]), line 334). The same atom therefore gets the same velocity whatever the number of processors, but a change of units changes its velocity (page A.4).

if (momentum_flag) zero_momentum(); step 2 if (rotation_flag) zero_rotation(); step 3 t = temperature->compute_scalar(); step 4: measure T rescale(t,t_desired); multiply v by √(t_desired/t)

2 · Predict by hand: the kinetic energy at step 0

LAMMPS defines the temperature from the kinetic energy and the number of degrees of freedom. For N atoms in 3D with the total momentum removed, there are 3N − 3 of them:

T = 2 KE / (dof × k_B) dof = 3N − 3 = 3 × 4000 − 3 = 11997 after step 4, T = 3 exactly, so KE = dof × T / 2 = 11997 × 3 / 2 = 17995.5 (total, in ε) the velocity of the centre of mass after step 2: 0 (to rounding, about 10⁻¹⁷)

With mom no, step 2 is skipped: the centre of mass keeps a small random velocity. No internal force can change it (forces between atoms cancel in pairs), so the prediction is that it stays exactly the same at every later step.

3 · Modify and observe

Choose the options, run 100 steps of melt, and compare the vx distribution at step 0 with the Maxwell bell curve at T = 3, and at step 100 with the Maxwell curve at the temperature LAMMPS prints for step 100.

dist mom rot seed

What the runs show

Every option gives KE = 17995.5 at step 0: step 4 always rescales to the exact target.
uniform starts flat, clearly not the Maxwell curve. After 100 steps (0.5 τ) collisions have already turned it into a bell curve. It is narrower than the T = 3 curve, because melting has turned kinetic energy into potential energy: the temperature has dropped to about 1.67, and the step-100 histogram follows the Maxwell curve for that temperature. The starting distribution is forgotten quickly; the temperature it sets is not.
mom no: the centre of mass moves at a small constant velocity, identical at steps 0 and 100 to every digit. The whole system drifts; with a thermostat this drift can grow into the "flying ice cube" of Lab 3.7. Keep the default.
A different seed: the same KE at step 0, different numbers at step 100. Runs with several seeds are the honest way to estimate error bars (Lab 6.8).

Your turn

1. What KE does LAMMPS print at step 0 for 500 atoms created at T = 1.5?

dof = 3 × 500 − 3 = 1497, KE = 1497 × 1.5 / 2 = 1122.75 (total, with thermo_modify norm no).

2. In a 2D run with 9 atoms (page A.3) at T = 0.72, what is the total KE at step 0?

3. Why are the random numbers divided by √m before rescaling?