Lab: two atoms by hand

The setup

Two argon atoms sit on a line. We use reduced units (Lesson 2.1): lengths in σ (3.4 Å for argon), energies in ε (the depth of the well), masses in units of the argon mass, so m = 1. One time unit is about 2.2 ps.

The atoms are slightly squeezed, so they will push apart and then vibrate like two balls on a spring. Nothing is added or taken away, so the total energy must stay constant. That is our test of whether a recipe can be trusted.

Step 1

The energy between two atoms

Why: MD needs a rule for how atoms push and pull. Real atoms repel strongly when squeezed and attract weakly when a little apart. The Lennard-Jones formula captures both with just two numbers (σ and ε).

In reduced units the formula is $U(r) = 4\left(r^{-12} - r^{-6}\right)$. Work out r⁻⁶ first, then square it to get r⁻¹²:

What the numbers say: The energy is lowest at r = 1.1225 (exactly −1 ε), rises steeply below r = 1 and fades to almost zero by r = 2.5. That last fact is why we can ignore far-away pairs (cutoffs, Lesson 2.5). Limitation: one formula for simple atoms only; molecules and ions need more terms (Chapters 4 and 5).

Step 2

From energy to force

Why: Newton's law needs forces, not energies. A force always pushes "downhill": it is minus the slope of the energy. Getting the force wrong is the most common bug in MD codes, so we also check it.

The slope of rⁿ is n·rⁿ⁻¹. Applying this to both terms:

Check with the slope measured directly between two nearby points (a "finite difference"):

What the numbers say: The force is large and repulsive when the atoms are squeezed, zero at the minimum, and weakly attractive beyond it. The finite-difference check costs one line and catches sign errors and missing factors.

Step 3

Newton's law for the pair

Why: To move atoms we need their accelerations, a = F/m. The two atoms feel equal and opposite forces (Newton's third law), so their total momentum stays zero.

What the numbers say: In every table below, v1 = −v2 exactly. If your own code ever gives a drifting total momentum, the forces are not equal and opposite (a bug).

Step 4

Recipe 1: Euler, three steps by hand

Why: The simplest idea: move each atom with its current velocity, and change the velocity with its current acceleration, both for the whole step. Let us see what happens to the energy.

What the numbers say: The energy rises by about 0.007 every step. Euler uses the force from the start of the step for the whole step, but the force is getting smaller as the atoms separate, so every kick is a little too big. The error always has the same sign, so it piles up.

Step 5

Recipe 2: velocity Verlet, three steps by hand

Why: Verlet uses the average of the old and new forces for the velocity, like the trapezoid rule for areas. It needs the same single force calculation per step as Euler, so it costs the same.

What the numbers say: After three steps the Verlet energy error is 0.0005, against 0.022 for Euler: about 40 times smaller for the same work. The small errors also go up and down instead of piling up (next step).

Step 6

Many steps: drift or wobble?

Why: Three steps do not prove anything about a whole simulation. We let the computer repeat our hand recipe for about three vibrations (2 time units) and compare different time steps.

What the numbers say: Euler adds energy steadily until the atoms fly apart: an unphysical result from a numerical error. Verlet keeps the energy within 0.2 % for dt = 0.01, forever. Its limitation: the energy is not exact (it wobbles), and the wobble grows as dt². This is why every MD code uses Verlet or its twin, leap-frog.

Step 7

How large may the time step be?

Why: A larger dt means fewer steps, so a cheaper simulation. But the time step must resolve the fastest motion in the system. We estimate that motion by hand and then test it.

Near the minimum the LJ well looks like a spring. Its stiffness is the curvature of U there, and two atoms of mass 1 vibrate like one particle with half the mass (the reduced mass μ = 1/2):

What the numbers say: For argon, dt = 0.01 τ ≈ 22 fs is fine. A water O–H bond vibrates with a period of about 9 fs, so water needs dt ≈ 0.5–1 fs, or constraints that freeze the bond (Lesson 4.4). The time step is set by the fastest motion, not by what you want to study.

Summary

The same steps with forces from many atoms, in three dimensions, are a complete MD program. Everything later in the course adds to this loop; nothing replaces it.