lammps-workflow

LAMMPS Workflow

Description

This skill covers end-to-end molecular dynamics simulations using LAMMPS (Large-scale Atomic/Massively Parallel Simulator): input script anatomy, force-field selection and deployment, ensemble control, equilibration protocol, production runs, restart logic, and parallel execution. LAMMPS is the dominant open-source MD code for materials science applications and the standard deployment target for ML interatomic potentials. Invoke this skill when setting up, running, diagnosing, or extending a LAMMPS simulation of any material class from metals to polymers to reactive systems.

Domain Context

LAMMPS integrates Newton's equations of motion for a set of interacting particles under a chosen force field. The equations are solved numerically using a symplectic integrator (Velocity Verlet by default), which conserves a shadow Hamiltonian close to the true Hamiltonian and preserves time-reversibility. Key physical approximations:

• Force field validity domain: Every pair_style encodes specific physical assumptions. EAM is validated for FCC metals near equilibrium; it fails for surfaces, defects with large distortions, or alloys outside the training composition range. Tersoff is validated for covalent semiconductors; it fails for amorphous structures far from the fitting database. ReaxFF can model bond breaking/forming but has known failures for unusual coordination environments. The force field, not the integrator, is the dominant source of error.
• Classical vs. quantum nuclei: LAMMPS integrates classical equations of motion. Zero-point energy and nuclear quantum effects are neglected. This is a good approximation for heavy atoms (≥ C) above ~200 K, but fails for light atoms (H, He, Li) at low temperature and for tunneling-dominated processes. [EXPERT REVIEW NEEDED for precise crossover conditions]
• Boundary conditions: LAMMPS supports periodic (p), shrink-wrapped (s), and fixed (f) boundaries. Periodic boundaries impose artificial periodicity; the simulation box must be large enough that a particle does not interact with its own image through the cutoff.
• Timestep and numerical stability: The timestep must resolve the fastest vibrational frequency in the system. A rule of thumb is dt ≤ period / 20. For stiff bonds (O-H, C-H) in molecular force fields, this requires dt ≤ 0.5 fs or constrained bonds (SHAKE/RATTLE). For metals with EAM, 1–2 fs is typically safe. For ReaxFF with short-range repulsion, 0.1–0.5 fs is required. Exceeding the stability limit produces exponentially growing kinetic energy — the "energy explosion" failure mode.
• Ergodicity and sampling: A single MD trajectory samples the microcanonical or canonical ensemble only if it is ergodic on the timescale of the run. Rare events (phase transitions, defect migration, folding) may not be sampled. Enhanced sampling methods are required for free energy calculations.