LAMMPS from scratch · 25 min

Run LAMMPS, change it, check it by hand

LAMMPS is the program most researchers use for molecular dynamics. Here it runs inside your browser. You run its classic first example, change one line at a time, and get every number LAMMPS prints at step 0 with a calculator, to the last digit.

What LAMMPS is, and how an input script tells it what to do
Run the official melt example live and read its output
Predict the effect of four changes by hand, then check them against LAMMPS

What LAMMPS is

LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) is a free molecular dynamics program from Sandia National Laboratories. Inside, it does what you built in Chapters 1 to 3: it computes forces, moves atoms with velocity Verlet and controls temperature. It does this for anything from 4 atoms to billions, on one laptop or thousands of processors.

You do not write a program to use LAMMPS. You write an input script: a plain text list of commands. LAMMPS reads it from top to bottom, carries out each command, and prints a log with the results.

The LAMMPS on this page is the real program, compiled to WebAssembly so that it runs inside your browser. Nothing is sent to a server. Its version is printed in the log when it runs.

How every page in this section works

1

Understand

Read the input or the code, one line at a time.

2

Modify

Change one thing, and predict the new result by hand first.

3

Observe

Run LAMMPS and see what changed.

4

Match

Compare LAMMPS with your hand calculation, digit by digit.

1 · Understand: the melt example, line by line

This script ships with LAMMPS (examples/melt/in.melt). It builds a Lennard-Jones crystal and gives it so much kinetic energy that it melts.

Line	What it does	Number you can check
`units lj`	Lennard-Jones units: σ = ε = m = 1 (Lesson 2.1)	energies in ε, lengths in σ
`atom_style atomic`	atoms are points: position, velocity, type
`lattice fcc 0.8442`	an FCC crystal with 0.8442 atoms per σ³	cube edge a = (4/0.8442)^1/3 = 1.6796
`region box block 0 10 0 10 0 10`	a box 10 cells long in each direction	L = 10a = 16.796
`create_box 1 box`	one atom type, in that box
`create_atoms 1 box`	fill the box: 4 atoms per FCC cell	N = 4 × 10³ = 4000
`mass 1 1.0`	mass of type 1	m = 1
`velocity all create 3.0 87287 loop geom`	random velocities, scaled to exactly T = 3.0; 87287 is the random seed	T = 3 at step 0
`pair_style lj/cut 2.5`	Lennard-Jones force, cut at 2.5σ (Lesson 2.5)
`pair_coeff 1 1 1.0 1.0 2.5`	ε = 1, σ = 1, cutoff 2.5 for type 1 with type 1
`neighbor 0.3 bin`	neighbour lists reach 2.5 + 0.3 = 2.8 (Lesson 2.6)
`neigh_modify every 20 delay 0 check no`	rebuild the lists every 20 steps	250/20 → 12 rebuilds
`fix 1 all nve`	move all atoms with velocity Verlet (Lesson 1.5)
`thermo 50`	print a line of results every 50 steps	steps 0, 50, …, 250
`run 250`	do 250 steps

2 · Run it in your browser

Press Run. The first time, your browser downloads LAMMPS (about 3 MB). You can edit the script before running.

Full log

3 · Match: step 0 by hand

At step 0 every atom still sits on the perfect crystal, so every atom has exactly the same neighbours. Within the cutoff 2.5 they form four shells, at distances r = a√(k/2) for k = 1, 2, 3, 4.

shell k	atoms n	distance r	r⁻⁶	energy per atom ½·n·4(r⁻¹² − r⁻⁶)
1	12	1.18765	0.356337	−5.5046614
2	6	1.67960	0.044542	−0.5106972
3	24	2.05708	0.013198	−0.6251271
4	12	2.37531	0.005568	−0.1328823
5	24	2.65569		beyond 2.5: not counted
sum	54			−6.7733681

A shortcut makes this a two-line calculation. Because a³ = 4/ρ, every r⁻⁶ is a simple fraction of one number, x = ρ²/2:

r_k⁶ = (a² k/2)³ = 2k³/ρ² so r_k⁻⁶ = x/k³, r_k⁻¹² = x²/k⁶ x = 0.8442²/2 = 0.35633682 x² = 0.12697593 S3 = 12/1 + 6/8 + 24/27 + 12/64 = 13.8263889 S6 = 12/1 + 6/64 + 24/729 + 12/4096 = 12.1296015 E_pair = ½ Σ n·4(r⁻¹² − r⁻⁶) = 2(x² S6 − x S3) = 2(1.54016742 − 4.92685145) = −6.7733681 KE per atom = (3/2) T (N−1)/N (3 degrees of freedom of = 1.5 × 3 × 3999/4000 the whole crystal's motion = 4.4988750 are removed) TotEng = −6.77336805 + 4.49887500 = −2.2744931 Press = ρT(N−1)/N + (ρ/3)·½ Σ n·r·F(r), r·F = 24(2r⁻¹² − r⁻⁶) = 0.8442 × 3 × 0.99975 + 4ρ(2x² S6 − x S3) = 2.53196685 + 3.3768 × (3.08033484 − 4.92685145) = 2.53196685 − 6.23531727 = −3.7033504

Run LAMMPS above to fill in its numbers.

4 · Modify, predict, observe, match

Each change below edits one or two lines. Do the prediction by hand first (it is worked out for you), then press the button: it changes the script, runs LAMMPS for step 0 only, and compares.

A. Cut the force shorter: 2.5 → 2.2

pair_style lj/cut 2.2 pair_coeff 1 1 1.0 1.0 2.2

Predict by hand:

the 4th shell (r = 2.3753) is now beyond 2.2: drop it S3 = 12 + 6/8 + 24/27 = 13.6388889 S6 = 12 + 6/64 + 24/729 = 12.1266718 E_pair = 2(0.12697593 × 12.1266718 − 0.35633682 × 13.6388889) = 2(1.53979542 − 4.86003829) = −6.6404857 TotEng = −6.64048574 + 4.49887500 = −2.1416107 Press = 2.53196685 + 3.3768 × (3.07959084 − 4.86003829) = 2.53196685 − 6.01221495 = −3.4802481

B. Weaker attraction: ε = 1.0 → 0.5

pair_coeff 1 1 0.5 1.0 2.5

Predict by hand:

every energy and every force is multiplied by ε E_pair = 0.5 × (−6.77336805) = −3.386684 TotEng = −3.38668403 + 4.49887500 = 1.112191 Press = 2.53196685 + 0.5 × (−6.23531727) = 2.53196685 − 3.11765864 = −0.58569179 (the kinetic parts do not depend on ε)

C. Squeeze the crystal: ρ = 0.8442 → 0.9

lattice fcc 0.9

Predict by hand:

x = ρ²/2 = 0.81/2 = 0.405 x² = 0.164025 (the 5th shell, 1.6441 × √(5/2) = 2.5996, is still outside 2.5) E_pair = 2(0.164025 × 12.1296015 − 0.405 × 13.8263889) = 2(1.98955789 − 5.59968750) = −7.2202592 TotEng = −7.22025922 + 4.49887500 = −2.7213842 Press = 0.9 × 3 × 0.99975 + 3.6 × (3.97911577 − 5.59968750) = 2.69932500 − 5.83405822 = −3.1347332

D. Start colder: T = 3.0 → 1.0

velocity all create 1.0 87287 loop geom

Predict by hand:

step 0: atoms still on the lattice, only velocities changed E_pair = −6.7733681 (does not depend on velocities) KE/N = 1.5 × 1.0 × 0.99975 = 1.49962500 TotEng = −6.77336805 + 1.49962500 = −5.2737431 Press = 0.8442 × 1.0 × 0.99975 + (−6.23531727) = 0.84398895 − 6.23531727 = −5.3913283

Your own change

Edit the script in step 2 yourself (for example lattice fcc 0.8 or a cutoff of 1.5) and press Run: the page works out the hand prediction for your script and compares it with LAMMPS.

What the four changes teach

Change	What moved	Why
A. cutoff 2.2	E_pair rises by exactly 0.1328823	the 12 atoms of the 4th shell no longer count: the reported energy depends on the cutoff, so papers must state it
B. ε = 0.5	E_pair halves, pressure changes, kinetic parts do not	ε multiplies every pair energy and force, never the velocities
C. ρ = 0.9	every shell moves closer; only x changes	all distances scale with a = (4/ρ)^1/3
D. T = 1.0	TotEng and Press move; E_pair does not	at step 0 the positions are the same; only velocities changed

Check every number in Python

Coming next: changing the C++ code itself

Here you changed the input. Later pages change the source code of LAMMPS with the same four steps: read the C++ lines that compute the force, change them (for example the 12 in r⁻¹² to a 9), predict the new energy by hand, and compare.

A browser cannot compile LAMMPS itself (compiling it took 4½ minutes on our 2-processor computer, with a compiler far too large to download into a web page). So each code-change page will come with the modified LAMMPS already compiled: you will see the exact lines changed, run the original and the modified program side by side here, and match the difference by hand. You can also rebuild it on your own computer with the steps in Part 0.

Your turn

1. Change the density to 0.8. Is the 5th shell inside the cutoff now? Predict E_pair.

2. Set the cutoff to 1.5 in both lines. Which shells remain, and what is E_pair?

3. In change B the total energy became positive. Does that mean the crystal flies apart?

The theory behind this page

Lesson 7.1 · From toy code to LAMMPS and GROMACS