Lab: Monte Carlo vs MD by hand

The setup

The Metropolis recipe for one trial move:

No velocities, no forces, no time step. Systems: 3 atoms on a line, a harmonic well, the 9-atom 2D LJ box of earlier labs (k_BT = 0.72), and the double well of Lab 7.4.

Step 1

Metropolis moves by hand

Why: Moves that lower the energy are always accepted; moves that raise it are accepted sometimes, with exactly the Boltzmann factor. Three LJ atoms on a line are enough to see the rule at work.

What the numbers say: A small uphill move is accepted 92 % of the time, a bigger one 57 %; a move into a much higher energy would almost never be accepted. Counting the old state again after a rejection is essential: it is what gives low-energy states their larger weight. Leaving it out (only recording accepted states) is a classic bug that biases every average.

Step 2

Does Metropolis really give Boltzmann?

Why: The acceptance rule is built so that, in the long run, each state is visited in proportion to e^(−U/k_BT) (detailed balance: the flow from A to B equals the flow from B to A). Test it on a harmonic well, where every average is known.

What the numbers say: Even the tails (|x| > 2, visited 4.6 % of the time) come out right. MC needs only energies, never forces, which makes it the natural method for potentials without smooth derivatives (hard spheres, lattice models, discontinuous cutoffs) and for moves that have no physical path at all (step 6).

Step 3

The same averages from MC and MD

Why: MC and MD sample the same canonical ensemble, so equilibrium averages must agree within error bars. This is a strong check of both codes, and it brings out details that are easy to miss.

What the numbers say: MC and MD agree within their error bars. Two details matter. First, MD with truncated forces samples the energy that is shifted to zero at r_c; MC must use the same shifted energy to be the same model (here the difference is within the errors, but for a short cutoff it is not). Second, our first MD run disagreed badly (−1.855, P = 0.31): the initial velocities had a net momentum, the thermostat kept part of the kinetic energy in the drift of the whole box, and the atoms were effectively colder. Comparing with MC exposed the bug. MC gives the kinetic part of the pressure exactly, so its pressure is less noisy (± 0.014 vs ± 0.046).

Step 4

Step size and efficiency

Why: The only MC parameter is the maximum step δ. Tiny steps are almost always accepted but change little; huge steps change a lot but are almost always rejected. The right measure is not the acceptance but how many sweeps it takes to get an independent sample, the statistical inefficiency g of Lab 6.8.

What the numbers say: Efficiency is flat and best for acceptances between about 25 % and 60 % (g ≈ 10 sweeps), and drops by a factor of 15 for timid steps (94 % acceptance) and by 3 for reckless ones (8 %). The old rule "aim for 50 %" is a rough guide; measuring g is the real test. For this small atomic system MC even decorrelates the energy with less work than MD (10 sweeps vs 200 steps). For molecules with bonds and for large systems, single-atom moves must be tiny to be accepted, and MD, which moves all atoms together along the forces, usually wins; that is why liquids are mostly simulated with MD.

Step 5

What MC cannot do: dynamics

Why: An MC sweep is not a unit of time. Its "trajectory" moves however the trial moves let it. Compute a mean squared displacement per sweep anyway and see what it means.

What the numbers say: The MC "diffusion coefficient" changes tenfold with the step size, a parameter we chose for efficiency: it says nothing about the liquid. MC gives equilibrium averages (energies, pressures, g(r), free energies) but no diffusion, viscosity, conductivity, lifetimes or spectra. Those need real dynamics, i.e. MD (or kinetic Monte Carlo with physical rates, a different method). The MD value, 0.0599, matches Lab 6.7.

Step 6

What MD cannot do: jump

Why: MD must follow a continuous path, so a barrier of 5 k_BT is crossed rarely. MC can propose any move whose acceptance it can compute, including a jump straight to the other well. In the double well U(s) = 5(s² − 1)², a "flip" s → −s has ΔU = 0.

What the numbers say: MD crossed once in 400 000 steps, so its 0.532 is luck. MC with flip moves crosses tens of thousands of times and nails 0.501. This is where MC shines: moves with no physical path (swapping two molecules, inserting and deleting particles in the grand canonical ensemble, changing the box volume, configurational-bias moves for chains, cluster moves near critical points). The limitation: designing such moves needs insight into the system, and they give no dynamics. Hybrid schemes combine both (MD for dynamics-like sampling plus MC swaps or insertions; hybrid Monte Carlo uses short MD trajectories as trial moves).

Summary