Point a computer at the night sky and you can grow your own galaxy. Give it a hundred thousand stars, tell each one that it is pulled by the gravity of every other, and let time crawl forward in tiny steps. Spiral arms curl out of nothing; clusters form, merge, and fling stray stars into the dark. Point that same machinery at a droplet of argon instead — a few thousand atoms nudging one another — and you watch a liquid boil, a crystal melt, a protein fold. It is the same program. Only the force law and the length scale change.

That program has a name in two dialects. Astronomers call it an N-body simulation; chemists and materials scientists call it molecular dynamics (MD). Both do exactly one thing: take many interacting particles and push them all forward together by integrating Newton's laws. This page is about how to do that well — because with millions of steps stacked end to end, a mediocre integrator quietly destroys your world, and a good one keeps it alive for as long as you care to watch.

The N-body problem: everybody pulls on everybody

Take N particles with positions \mathbf{r}_1, \dots, \mathbf{r}_N. The rule of the game is simple to state: each particle feels the sum of the forces from all the others. For gravity, the pull on particle i from particle j is

which points from i toward j and falls off as 1/r^2 in strength (the cube in the denominator is one power for the 1/r^2 law and one to normalise the direction vector). The total force on i, and hence its acceleration \mathbf{a}_i = \mathbf{F}_i/m_i, is the sum over every other particle:

Here is the sting in the tail. To take a single step you must, for every one of the N particles, add up contributions from the other N-1. That is N(N-1)/2 distinct pairwise interactions per step — the cost grows as N^2. Double the number of stars and you quadruple the work; a thousand-fold bigger simulation is a million times slower. This O(N^2) wall is the central obstacle of the whole field, and taming it (see the tree-code digression below) is what makes billion-body cosmology possible.

For N = 2 Newton solved this by hand and got ellipses. For N = 3 there is no general closed-form solution — the three-body problem is famously chaotic. So from three particles upward, the computer is not a convenience; it is the only way through.

The heart of the matter: velocity-Verlet

Once you can compute every particle's acceleration, you need to march the positions forward in time. The obvious choice — plain Euler, "new position = old position + velocity × step" — is a disaster over long runs. It has a fatal habit: it steadily leaks or pumps energy. An Euler-integrated orbit spirals slowly outward (or inward) not because the physics says so, but because the method quietly adds a little false energy every step. Run it for a million steps and your tidy ellipse has unwound into a spiral. For a simulation whose whole point is to run for a very long time, that is fatal.

The workhorse of MD and N-body work is instead velocity-Verlet. It updates position and velocity in a carefully paired dance, using the acceleration at both the old and the new position:

Why is this modest-looking scheme the near-universal choice, beating the fancier Runge–Kutta methods that dominate elsewhere? Four reasons, and they matter:

• It is symplectic. Verlet preserves the geometric structure of the equations of motion, so its energy error stays bounded — it wobbles around the true value forever instead of drifting away. Your orbit stays an orbit for a billion steps.
• It is time-reversible. Run it forward then flip the velocities and run back, and you return exactly to the start (to round-off). The physics has that symmetry; so does the method.
• It is cheap. Just one force evaluation per step. Since computing all those O(N^2) pairwise forces is where nearly all the time goes, a method that asks for the force four times a step (like classical RK4) is four times as slow for the same \Delta t.
• It is second-order accurate. Halve the step and the per-step error drops by a factor of four — accurate enough, and paired with the bounded-energy property, it is the sweet spot for long simulations.

The tension is worth stating plainly: RK4 is more accurate on any single step than Verlet. But accuracy per step is the wrong thing to optimise for a run of millions of steps. What you want is for the total energy not to drift — and there, cheap symplectic Verlet crushes expensive non-symplectic RK4. The next block lets you watch the difference with your own eyes.

Verlet vs. Euler: watch the energy

Here is a body orbiting a fixed mass (units where GM = 1), launched on a circular orbit of radius 1. The total energy per unit mass, E = \tfrac{1}{2}v^2 - 1/r, should stay pinned at -0.5 forever. We integrate the same orbit two ways and print the energy as it goes. Watch Verlet hold steady while Euler slides away.

The Verlet energy jitters in the fifth decimal place and comes right back; the Euler energy marches off in one direction and never returns. Multiply that drift by the millions of steps a real simulation takes and the Euler orbit is long gone. This single output is the whole argument for symplectic integrators, made concrete.

From two bodies to many

Scaling up to a genuine N-body system changes nothing about the integrator — it changes only the force routine, which now has to loop over all pairs. Here are three equal masses released near one another (with a little "softening" added to the distance so a close encounter can't blow up to infinity), advanced with the very same velocity-Verlet dance. Notice the inner double loop: that is the O(N^2) cost in the flesh, doing exactly N(N-1)/2 = 3 pairwise force calculations per step.

The three bodies swing around their common centre in a genuinely three-body dance — no formula could have told you these coordinates in advance. To simulate a galaxy you change three numbers: N becomes a hundred thousand, the initial positions become a disc, and the loop becomes the bottleneck. Everything else is identical.

Molecular dynamics: atoms and the Lennard–Jones potential

Swap gravity for the force between two neutral atoms and you have molecular dynamics. Atoms do not pull on each other like planets. Push two of them together and they repel ferociously (their electron clouds refuse to overlap); pull them apart and they feel a gentle attraction (the fleeting van der Waals stickiness). The classic model that captures both in one tidy formula is the Lennard–Jones potential:

The (\sigma/r)^{12} term is the steep repulsive core — the wall an atom hides behind; the -(\sigma/r)^{6} term is the softer, longer-ranged attractive tail. Here \varepsilon sets the depth of the well (how sticky) and \sigma sets the distance scale (how big). The curve dips to a minimum — the comfortable spacing atoms settle into — at r = 2^{1/6}\sigma \approx 1.122\,\sigma, where the potential equals exactly -\varepsilon. The force along the line joining the atoms is just the downhill slope of this curve:

In a real MD run you drop a few thousand of these atoms into a box, give them random starting velocities, and integrate — with velocity-Verlet, of course. Two ideas make such a run physical:

• Temperature is average kinetic energy. The temperature of the sample is not an extra ingredient you impose; it simply is the mean kinetic energy of the atoms, \tfrac{3}{2}k_B T = \langle \tfrac{1}{2}mv^2 \rangle. Heat the box and the atoms jiggle harder; cool it and they settle into a lattice. "Melting" is nothing more than the jiggling getting violent enough to break the Lennard–Jones cages.
• Periodic boundary conditions. A few thousand atoms is a speck; to imitate bulk matter without walls, an atom that drifts out the right face of the box reappears at the left, and it feels forces from the nearest copies of its neighbours. The little box thus pretends to be an infinite crystal, dodging the artefacts a hard container would create.

An atom pair, oscillating in its well

Let us make that Lennard–Jones curve move. Take a single pair of atoms in one dimension (reduced units \varepsilon = \sigma = 1), start them stretched apart at r = 1.5 and release from rest. They should fall inward toward the bottom of the well near 1.122, overshoot, and swing back — a vibrating molecule. The energy sloshes between kinetic and potential, but the total holds steady, because Verlet is on the job.

Follow the r column: it dives from 1.5 down through the minimum and back out — an oscillation, exactly like a mass on a spring, only the "spring" is the curved floor of the Lennard–Jones well. The \text{KE} and \text{PE} columns trade places while the E column barely twitches. A whole vibrational spectrum of matter is hiding in that one little loop.