Lesson 7.1 · 15 min

From toy code to LAMMPS and GROMACS

Everything in your toy code has a one-line counterpart in LAMMPS and GROMACS: the lattice, the pair potential, the neighbour list, the integrator and the thermostat. This lesson runs the same Lennard-Jones liquid in our NumPy code and in LAMMPS and gets identical numbers, then shows how to translate between unit systems and how a real simulation workflow is organized.

How each part of a toy MD code maps onto LAMMPS and GROMACS input commands
Unit systems (LAMMPS lj, real, metal; GROMACS nm, ps, kJ/mol) and how to convert between them
The standard workflow (build, minimize, NVT, NPT, production, analysis) and basic sanity checks

The idea

Research codes such as LAMMPS and GROMACS do exactly what our toy code does, only much faster and with many more options. Every piece maps onto an input command:

Starting structure: lattice + create_atoms or read_data (LAMMPS); a .gro/.pdb file (GROMACS).
Force field: pair_style lj/cut 2.5 + pair_coeff (LAMMPS); [ atomtypes ] in a .top file plus cutoffs in the .mdp (GROMACS).
Neighbour list (Lesson 2.6): neighbor 0.3 bin; GROMACS sets its Verlet buffer automatically.
Integrator: fix nve (velocity Verlet); GROMACS uses leap-frog (integrator = md), equivalent for most purposes.
Thermostat and barostat (Chapter 3): fix nvt/fix npt (Nosé–Hoover); tcoupl = v-rescale, pcoupl = C-rescale.

Same model, same numbers. If you give LAMMPS exactly the same positions, velocities, potential, cutoff and time step as our NumPy code, you get the same energies to 10 significant figures (code tab, Part 2). That is a powerful way to test your own code, or to understand what a command really does.

Units are the classic trap. LAMMPS has several unit systems: lj (reduced: σ = ε = m = 1), real (Å, fs, kcal/mol, atm) and metal (Å, ps, eV, bar), among others. GROMACS always uses nm, ps, kJ/mol and bar. For argon (σ = 3.405 Å, ε/k_B = 119.8 K), τ = 2.156 ps, so Δt* = 0.005 is 10.8 fs, T* = 0.722 is 86.5 K and ρ* = 0.8442 is 1.419 g/cm³ (Part 3).

A real workflow has stages: build the system (Packmol, gmx insert-molecules, moltemplate, fftool), minimize the energy, equilibrate in NVT, then NPT to find the density, then a production run, then analysis (Chapter 6). At each stage check the basics: energy conservation in NVE, temperature and pressure fluctuating around their targets, a sensible density, and a sensible g(r).

Speed. Real codes split the box into domains, one per CPU core (domain decomposition), use neighbour lists and PME (Lesson 5.3), and run on GPUs. That is how they handle millions of atoms where our NumPy code handles hundreds.

Key idea. A research code is your toy code with better engineering. Map each concept to its command, keep the units straight, start from a working example, and validate against something you understand.

The algorithm

1. build      : coordinates + topology/force field   (Packmol, gmx insert-molecules, moltemplate)
2. minimize   : LAMMPS minimize / GROMACS integrator = steep   -> remove overlaps
3. NVT        : heat to the target T (short, e.g. 100 ps)
4. NPT        : equilibrate density (check that the box volume has levelled off)
5. production : NVT or NPT, save coordinates (and velocities) at the interval your analysis needs
6. analysis   : g(r), MSD, conductivity ... (Chapter 6) with error bars from blocks or replicas
sanity checks : NVE drift, , , density vs experiment, g(r) shape, no atoms lost

LAMMPS: one input script (in.*) runs everything; lmp -in in.lj or mpirun -np 8 lmp -in in.lj. GPU and KOKKOS packages accelerate it.
GROMACS: gmx grompp combines .mdp + .gro + .top into a .tpr; gmx mdrun runs it, using GPUs automatically when available.
Keep the input files, code version and random seeds with every data set, so that a run can be reproduced.

See it inside LAMMPS

The same idea in the real code, run live in your browser and checked by hand.