LAMMPS from scratch · 25 min

Forces and neighbours in the script

Five commands control how forces are computed: pair_style, pair_coeff, pair_modify, neighbor and neigh_modify. You predict by hand what shifting the energy, adding a tail correction and mixing two atom types do, run each live, and then examine how the melt example's neighbour-list settings behave in this particular run.

What pair_style, pair_coeff and pair_modify (shift, tail, mix) change, predicted and checked
How many neighbours each atom has in the list, counted by hand
How neigh_modify trades speed against missed interactions ("dangerous builds")

1 · Understand: five commands

Command	Sets	In melt
`pair_style lj/cut 2.5`	the kind of force, and the default cutoff	Lennard-Jones, cut at 2.5
`pair_coeff 1 1 1.0 1.0 2.5`	ε, σ (and optionally the cutoff) for one pair of atom types	types 1–1: ε = 1, σ = 1
`pair_modify …`	`shift yes`: energy 0 at the cutoff; `tail yes`: add the part beyond the cutoff; `mix`: rule for unlike types	not used
`neighbor 0.3 bin`	the skin: neighbour lists reach cutoff + 0.3 = 2.8; bin: how they are built (Lesson 2.6)	2.8
`neigh_modify every 20 delay 0 check no`	when to rebuild the lists	every 20 steps, no checking

2 · Modify, predict, observe, match

Shift the energy to zero at the cutoff

pair_modify shift yes

Predict by hand:

without the shift, u jumps from u(2.5) to 0 at the cutoff: u(2.5) = 4(2.5⁻¹² − 2.5⁻⁶) = 4(0.0000168 − 0.004096) = −0.01631689 shift yes subtracts u(2.5) from every pair inside the cutoff: each atom has 54 neighbours within 2.5, each pair counted half: ΔE_pair = −½ × 54 × (−0.01631689) = +0.44055606 E_pair = −6.77336805 + 0.44055606 = −6.332812 TotEng = −6.33281199 + 4.49887500 = −1.833937 Press = −3.7033504 unchanged: a constant shift does not change the force

Add the missing tail beyond the cutoff

pair_modify tail yes

Predict by hand:

beyond 2.5 the liquid is assumed uniform (g = 1): E_tail = (8/3)πρ [ (1/3) r_c⁻⁹ − r_c⁻³ ] = 7.0723534 × (0.0000874 − 0.064) = −0.45201262 P_tail = (16/3)πρ² [ (2/3) r_c⁻⁹ − r_c⁻³ ] = 11.940961 × (0.0001748 − 0.064) = −0.76213470 E_pair = −6.77336805 − 0.45201262 = −7.2253807 TotEng = −2.27449305 − 0.45201262 = −2.7265057 Press = −3.70335042 − 0.76213470 = −4.4654851

Two atom types: mixing rules

Give type 2 its own ε = 0.5 and σ = 1.2 and say nothing about the pair 1–2. LAMMPS invents it with a mixing rule; lj/cut uses the geometric rule unless told otherwise. Two atoms, type 1 at x = 0 and type 2 at x = 1.2:

geometric : ε₁₂ = √(1.0 × 0.5) = 0.70710678 σ₁₂ = √(1.0 × 1.2) = 1.0954451 arithmetic: ε₁₂ = √(1.0 × 0.5) = 0.70710678 σ₁₂ = (1.0 + 1.2)/2 = 1.1 E = 4ε₁₂[(σ₁₂/r)¹² − (σ₁₂/r)⁶] fx on atom 1 = −24ε₁₂[2(σ₁₂/r)¹² − (σ₁₂/r)⁶]/r geometric : (σ/r)⁶ = 0.57870370 E = −0.68958673 fx = −1.2882389 arithmetic: (σ/r)⁶ = 0.59329219 E = −0.68248975 fx = −1.5655211

Force fields are built for one rule: OPLS uses geometric, AMBER and CHARMM use arithmetic (Lorentz–Berthelot). Using the wrong one changes every unlike interaction.

3 · The neighbour list, counted by hand

The list holds every pair closer than cutoff + skin = 2.8. In the perfect crystal at step 0 the shells inside 2.8 are 12 + 6 + 24 + 12 + 24 = 78 atoms (the fifth, at 2.6557, is inside the skin but outside the cutoff). Each pair is stored once (Newton's third law), so:

neighbours per atom = 78 / 2 = 39 total in the list = 4000 × 39 = 156000

4 · Observe: are melt's neighbour settings safe?

The list can miss an interaction if, between rebuilds, an atom moves far enough that a pair previously outside the list cutoff (2.8) enters the force cutoff (2.5). LAMMPS flags a dangerous build when its displacement check detects this condition at a rebuild. With check no, that check is disabled, so LAMMPS cannot report whether scheduled rebuilds were dangerous. Whether this schedule misses interactions depends on the timestep, skin, and actual particle motion; the estimate below indicates why this setup deserves caution.

neigh_modify	idea	list builds	dangerous	TotEng at 250
`every 20 delay 0 check no`	melt as shipped
`every 20 delay 0 check yes`	same schedule, but count the dangerous builds
`every 10 delay 0 check yes`	look twice as often
`every 1 delay 0 check yes`	look every step, rebuild only when needed

In the recorded run, the melt settings rebuilt the list 12 times, but with check no LAMMPS did not check or report dangerous builds. With every 1 delay 0 check yes, it rebuilt 42 times (about every 6 steps) and reported no dangerous builds. The total energy at step 250 was −2.2812174 for the original settings and −2.2797308 for the checked settings, versus −2.2744931 at step 0. These results are specific to this input and LAMMPS build; use a checked neighbor schedule when validating a simulation.

Your turn

1. With pair_modify shift yes, why does the pressure not change at all?

2. With a skin of 0.5 instead of 0.3, how many neighbours per atom does the step-0 list hold?

3. The tail correction lowers the pressure by 0.76. Is that small?