LAMMPS
LAMMPS, standing for Large-scale Atomic/Molecular Massively Parallel Simulator, is an open-source classical molecular dynamics (MD) simulation software package focused on materials modeling.[1] Developed primarily by researchers at Sandia National Laboratories, it integrates Newton's equations of motion to model ensembles of particles—ranging from atoms to coarse-grained or continuum entities—in liquid, solid, or gaseous states across 2D or 3D systems.[2] LAMMPS excels in simulating complex materials like metals, semiconductors, polymers, biomolecules, and granular systems, scaling from a few to billions of particles on single processors or parallel platforms including supercomputers.[3] Its design emphasizes efficiency through spatial decomposition, neighbor lists, and optimizations for short-range interactions, while supporting extensions via modular packages.[2] The development of LAMMPS originated in the mid-1990s under a Cooperative Research and Development Agreement (CRADA) between U.S. Department of Energy laboratories—Sandia National Laboratories and Lawrence Livermore National Laboratory—and industry partners including Cray Research, Bristol-Myers Squibb, and DuPont, led by Steve Plimpton at Sandia.[4] Initial versions were released in Fortran 77 (1999) and Fortran 90 (2001), incorporating capabilities from related codes like ParaDyn and GRASP, before the current C++ codebase was publicly released as open source on September 1, 2004, under the GNU General Public License version 2.[4] Since then, LAMMPS has undergone continuous enhancements, including the addition of many-body potentials like ReaxFF in 2006 and requirements for C++11 compatibility by 2020, with the latest stable release on July 22, 2025.[4][5] Primarily funded by the U.S. Department of Energy, it received an R&D 100 Award in 2018 for its innovations in scalable simulation.[1] Key features of LAMMPS include support for diverse particle types—such as point particles, finite-size atoms, dipoles, and rigid bodies—and a broad array of interatomic potentials, encompassing pairwise (e.g., Lennard-Jones), many-body (e.g., embedded atom method), bond/angle/cosine styles for polymers, long-range electrostatics via Ewald or particle-particle particle-mesh methods, and even machine-learning potentials.[6] It offers flexible simulation controls through input scripts with variables, loops, and conditional logic; various thermodynamic ensembles (e.g., NVE, NVT, NPT); integrators like velocity-Verlet; and boundary conditions for periodic, fixed, or shrink-wrapped systems.[6] Parallel execution leverages the Message Passing Interface (MPI) with optional OpenMP threading, GPU acceleration (via CUDA, HIP, OpenCL, or SYCL), and Intel Xeon Phi support, enabling high performance on distributed-memory clusters.[2] Additional capabilities include atom creation/deletion, diagnostics, and output in formats like XYZ or LAMMPS dump files, with extensibility through a C library interface, Python wrappers, and plugins for custom force fields or models.[6]History and Development
Origins
Development of LAMMPS began in the mid-1990s under a Cooperative Research and Development Agreement (CRADA) between Sandia National Laboratories, Lawrence Livermore National Laboratory, Cray Research, Bristol-Myers Squibb, and DuPont.[4] The project aimed to produce a large-scale parallel classical molecular dynamics (MD) code capable of simulating materials on emerging parallel computing architectures.[4] The coding effort was led by Steve Plimpton at Sandia National Laboratories, with the primary goal of creating a flexible MD simulation tool that could model materials without requiring extensive custom modifications to existing codes.[4][7] At the time, many MD packages were limited in scalability, often taking hours on vector supercomputers like the Cray Y-MP for simulations of 10,000 to 100,000 atoms over picoseconds, or inefficient on SIMD machines due to issues with global data structures and short-range forces.[8] LAMMPS was designed to address these constraints by prioritizing parallel efficiency for systems with rapidly changing neighbors, with an initial focus on polymer and granular simulations.[8] Early implementations were written in Fortran 77, emphasizing computational efficiency for large-scale atomic and molecular systems on parallel platforms such as the nCUBE 2 and Intel iPSC/860, achieving up to 90% parallel efficiency for systems with millions of atoms.[7][8] These pre-2000 versions laid the foundation for extensible force computations and neighbor lists tailored to short-range interactions common in materials modeling.[8] Over time, the code evolved into an open-source C++ version to enhance portability and community contributions.[4]Key Milestones and Releases
LAMMPS development concluded its initial phase under a Cooperative Research and Development Agreement (CRADA) in 1999 with the release of LAMMPS 99, the final version written in Fortran 77.[4] This marked the end of proprietary restrictions from the CRADA between Sandia National Laboratories, Lawrence Livermore National Laboratory, and industrial partners including Cray Research, Bristol-Myers Squibb, and DuPont.[4] In 2001, LAMMPS transitioned to Fortran 90 with the release of LAMMPS 2001, introducing improved memory management to enhance scalability for larger simulations.[4] A major architectural shift occurred in 2004 when LAMMPS was rewritten in C++ and released as open-source software under the GNU General Public License version 2 (GPLv2), enabling broader community contributions.[4] This version incorporated parallel computing elements from related codes, including ParaDyn for dynamic remapping, Warp for short-range interactions, and GranFlow for granular simulations.[4] By 2006, LAMMPS integrated capabilities from Aidan Thompson's GRASP code, expanding support for many-body interatomic potentials such as the Stillinger-Weber, Tersoff, and ReaxFF models to better handle complex material behaviors.[4] In March 2020, the stable release introduced a requirement for C++11-compatible compilers, aligning with modern software standards to facilitate advanced features while maintaining backward compatibility for core functionality. Recent releases include the stable version on July 22, 2025 (with an update on November 13, 2025), followed by a feature release cycle beginning on September 10, 2025, which implemented fundamental changes to the build system (raising minimum CMake to 3.20 and restricting the traditional build system) and code structure (including support for C++17 and removal of certain packages).[9] This release also built on prior enhancements, such as the addition of the dispersion/d3 pair style for D3 dispersion corrections with DFT-based machine learning potentials in the July 2025 stable release.[9] Biennial LAMMPS workshops support ongoing development, with the Ninth LAMMPS Workshop and Symposium held August 12-14, 2025, in Albuquerque, New Mexico, featuring tutorials, technical talks, and community discussions.[10]Contributors and Community
LAMMPS was primarily developed by Steve Plimpton at Sandia National Laboratories, who created the initial codebase in the mid-1990s and continues to lead its maintenance.[1][11] A key contribution came from Aidan Thompson, also at Sandia, who integrated features from his GRASP molecular dynamics code into LAMMPS starting in 2006, enhancing its support for reactive force fields and other advanced models.[4][12] Development and maintenance of LAMMPS have been primarily funded by the United States Department of Energy (DOE) through programs such as CRADA, LDRD, ASCI, ECP, and others, enabling Sandia-led efforts to sustain and expand the software.[13] This funding supports core development activities, including performance optimizations and integration of new computational methods.[13] Since its release as open-source software in 2004 under the GNU General Public License, LAMMPS has fostered global contributions through its GitHub repository, where users worldwide submit code enhancements, bug fixes, and new features.[4][14] As of May 2025, the core codebase consists of approximately 750 source files totaling around 165,000 lines of C++ code, with over 430 total contributors and a core development team of 12 from five institutions across the US and Europe.[12] The LAMMPS community is governed through collaborative channels including the MatSci.org forum (succeeding the lammps-users mailing list), GitHub issue trackers for bug reporting and feature requests, and biennial workshops that facilitate feedback and training.[15][16][17] These mechanisms ensure ongoing engagement, with workshops held roughly every two years in Albuquerque, New Mexico, to discuss advancements and user needs.[16] Recent modernization efforts, detailed in a 2025 publication, have applied contemporary software engineering practices to LAMMPS as a research code, emphasizing modular design with 92 optional packages for extensibility and comprehensive testing suites.[18] These include automated unit tests using GoogleTest (achieving about 80% file coverage and 45% line coverage), regression tests via the examples directory, and static analysis tools like Coverity Scan and CodeQL, all integrated with continuous integration systems such as Jenkins and Azure GitHub runners.[18] This refactoring promotes maintainability, with transitions to modern C++ standards (C++11 since 2020 and C++17 planned for late 2025) and a CMake-based build system replacing legacy tools.[18]Core Architecture and Capabilities
Parallel Computing and Performance
LAMMPS utilizes a spatial-decomposition algorithm to enable efficient parallel execution, dividing the simulation domain into rectangular subdomains assigned to individual MPI processes. This approach distributes atoms and computations across processors, with each process responsible for calculating forces on atoms within its subdomain, including interactions with neighboring subdomains via ghost atoms that extend the subdomain boundaries by a cutoff distance. The method minimizes communication overhead by localizing most computations and using MPI point-to-point messages for irregular data exchanges, such as atom migrations during time steps. This design, introduced in the original LAMMPS implementation, supports both 2D and 3D simulations and scales effectively with increasing processor counts.[19] The code supports distributed-memory parallelism through the Message Passing Interface (MPI) standard, allowing simulations to run across thousands of processors on high-performance computing clusters. For shared-memory systems, LAMMPS incorporates OpenMP threading to parallelize loops within each MPI process, particularly for force computations and neighbor list builds in supported styles from the OPENMP package. Hybrid parallelism combines MPI across nodes with OpenMP threads per node, offering flexibility for multi-core architectures; for instance, configurations like 4 MPI tasks per node with 4 threads each can optimize performance on systems with limited interconnect bandwidth. This hybrid mode is invoked via command-line options or input script commands, with benchmarks showing 5-20% speedups from single-threaded OpenMP over serial styles in certain scenarios.[19][20] To leverage accelerators, LAMMPS includes packages for GPU computing via CUDA for NVIDIA hardware, OpenCL for various GPUs, HIP for AMD GPUs, and SYCL for Intel GPUs, as well as support for Intel Xeon Phi co-processors through the INTEL package. The KOKKOS package provides performance portability by abstracting hardware-specific kernels into a single C++ codebase, enabling execution on CPUs, GPUs, and other devices without code modifications; it supports backends like CUDA, HIP, OpenMP, and SYCL. Integrations as of the 2025 releases (stable version 22 July 2025, with updates through November 2025) have enhanced KOKKOS for heterogeneous computing environments, including new features like accelerated fixes for electron stopping and improved multi-device utilization on modern supercomputers.[21][22][23][24] LAMMPS demonstrates strong scalability on exascale supercomputers, achieving parallel efficiencies exceeding 90% for weak scaling (system size proportional to processors), with recent benchmarks (as of 2025) showing performance scaling to over 8,000 nodes on systems like Frontier. Load-balancing is addressed through the balance command, which dynamically adjusts subdomain sizes and shapes—using grid-based or recursive bisection methods—to equalize particle counts or computational weights across processors, mitigating imbalances from uneven densities or hybrid interaction costs. This is particularly useful for interfaces or clustered systems, with rebalancing triggered when load variance exceeds a user-defined threshold, such as 20%. Overall, these optimizations enable LAMMPS to perform efficiently on petascale and exascale machines, with representative benchmarks showing sustained flop rates in the petaflops range for billion-atom systems.[25][26][23]Supported Particle Systems and Models
LAMMPS supports a wide range of particle representations through its atom_style command, enabling simulations of diverse physical systems at the atomic, molecular, and mesoscale levels.[27] The core atomic style provides the minimal attributes for point particles, such as position, velocity, and force, suitable for modeling simple atomic systems like metals, semiconductors, and liquids where particles interact as point masses.[27] For charged systems, the charge style extends this by including charge attributes, allowing simulations of ionic materials or electrolytes.[27] Molecular systems are accommodated by styles like bond, which adds connectivity for bonded interactions in bead-spring polymer models, and molecular, which incorporates angles, dihedrals, and impropers for apolar or uncharged molecules such as hydrocarbons or polymers.[27] The full style combines these with charges, supporting all-atom representations of biomolecules, proteins, and complex organic molecules in aqueous environments.[27] Coarse-grained representations, such as the oxdna style for DNA and RNA strands or template for templated molecular systems, enable mesoscale modeling of larger biomolecular assemblies by reducing detail while capturing essential dynamics.[27] Extended particle representations handle non-spherical or rigid bodies, including the sphere style for finite-size spherical particles with radius, density, and rotational degrees of freedom, commonly used in granular flows or colloidal suspensions.[27] The ellipsoid style adds shape and orientation parameters for aspherical particles, while body allows arbitrary rigid or deformable bodies with internal degrees of freedom, such as polyhedral objects or integrated finite element models.[27] Specialized 2D and 3D rigid body styles, like line for line segments and tri for triangulated surfaces, support simulations of meshed or polygonal particles in lower-dimensional or surface-based systems.[27] Granular and hybrid systems are facilitated by granular variants of the sphere style, incorporating friction and cohesion for modeling powders, sands, or bonded particle media, and the hybrid style, which combines multiple sub-styles to simulate heterogeneous systems, such as metals interfaced with polymers or atomic regions coupled to coarse-grained domains.[27] Additional extended styles include dipole for point dipole particles in polarizable models, spin for magnetic spins in ferromagnetic materials, and dielectric for surface polarization effects.[27] Mesoscopic methods feature styles like dpd for dissipative particle dynamics in soft matter flows, sph for smoothed particle hydrodynamics in fluids, and peri for peridynamic continuum damage modeling, all operating in classical frameworks without quantum mechanical treatments.[27] LAMMPS operates in both 2D and 3D geometries, with particle interactions defined through pair, bond, and other potentials as detailed in subsequent sections.[27]Force Fields and Interatomic Potentials
LAMMPS supports a wide array of force fields and interatomic potentials to compute atomic forces in molecular dynamics simulations, ranging from simple pairwise interactions to complex many-body and reactive models. These potentials enable modeling of diverse materials, including metals, semiconductors, polymers, and biomolecules, by approximating the potential energy as a function of atomic positions. The implementation allows users to select appropriate styles via thepair_style command, with parameters specified through pair_coeff directives.[28]
Pairwise additive potentials in LAMMPS form the foundation for non-bonded interactions, treating forces between atom pairs independently. The Lennard-Jones (LJ) potential, a 12-6 form modeling van der Waals attractions and repulsions, is widely used for noble gases and simple fluids; it is implemented with cutoff options like lj/cut for efficiency.[28] The Buckingham potential, employing an exponential repulsion and inverse-power attraction, suits ionic solids and alkali halides, available as buck.[28] The Morse potential, with its exponential form for bond stretching, applies to diatomic molecules and surface interactions, via morse.[28] Coulombic interactions, representing electrostatic forces, are handled pairwise with short-range cutoffs (coul/cut) or long-range summation using Ewald techniques (coul/long with kspace_style ewald) or particle-particle particle-mesh (PPPM) methods (pppm), which accelerate computations for periodic systems by solving Poisson's equation in Fourier space.[28][29]
Many-body potentials extend beyond pairwise sums by incorporating environmental dependencies, capturing effects like metallic cohesion and covalent bonding. The Embedded Atom Method (EAM), implemented as eam/alloy or eam/fs, models metals by embedding an atom in the electron density of neighbors, with the energy including a pairwise term and an embedding function; it was developed for fcc metals like copper and nickel.[28] The Modified Embedded Atom Method (MEAM), via meam, adds angular terms to EAM for directional bonding in alloys and hcp structures, enhancing accuracy for systems like titanium.[28] The Tersoff potential, using tersoff, employs bond-order terms dependent on local coordination for covalent materials such as silicon and diamond, balancing two- and three-body interactions.[28] For reactive systems, ReaxFF (reaxff) simulates bond breaking and formation in hydrocarbons and organics through variable bond orders and charge equilibration, without predefined connectivity.[28]
Specialized potentials address coarse-grained or anisotropic systems. Dissipative Particle Dynamics (DPD), via dpd, models mesoscale hydrodynamics in soft matter like polymers by adding conservative, dissipative, and random forces to coarse-grained beads, preserving momentum.[28] The Gay-Berne potential (gayberne) extends LJ for ellipsoidal particles, capturing orientational dependencies in liquid crystals through anisotropic attractions.[28] Water models like TIP3P (tip3p) and SPC (spce) represent rigid three-site molecules with partial charges on hydrogens and oxygen, plus LJ on oxygen, for biomolecular simulations; TIP3P uses a smaller charge and geometry to match liquid properties.[28]
LAMMPS facilitates combining potentials through hybrid styles (pair_style hybrid), overlaying interactions like LJ with EAM for mixed systems, or using pair_style hybrid/overlay for additive contributions. Additionally, integration with the Open Knowledgebase of Interatomic Models (OpenKIM) via pair_style kim allows access to a repository of user-contributed, verified potentials, including classical and machine-learned models, by specifying KIM model IDs.[30][31]
Machine learning potentials, treated as advanced force field types in LAMMPS, include Atomic Cluster Expansion (ACE), Gaussian Approximation Potentials (GAP), and Neural Network Potentials (N2P2 or Behler-Parrinello), enabling high-fidelity predictions for complex materials without explicit functional forms.[6]
Simulation Methods
Boundary Conditions, Ensembles, and Constraints
LAMMPS supports a variety of boundary conditions to define the simulation box geometry and particle interactions at its edges, enabling simulations in both two-dimensional (2D) and three-dimensional (3D) spaces. In 3D simulations, boundary conditions are specified independently for each pair of opposing faces (x, y, z directions) using theboundary command. The primary types include periodic boundaries (p), where particles exiting one face re-enter through the opposite face, allowing seamless interactions across the box; fixed boundaries (f), which maintain constant face positions and delete atoms that cross them (unless lost atoms are ignored via thermo_modify lost ignore); and shrink-wrapped boundaries (s), which dynamically adjust face positions to encompass all atoms, with a variant (m) enforcing a minimum distance based on the initial box size to prevent collapse. For 2D simulations, the z-dimension defaults to periodic boundaries to simulate a slab geometry. LAMMPS also accommodates triclinic (non-orthogonal) boxes through tilt factors (xy, xz, yz), where periodic boundaries account for shear deformations, and shrink-wrapped boundaries apply to tilted faces, facilitating simulations of anisotropic systems like sheared fluids.[32][33]
Thermodynamic ensembles in LAMMPS control the simulation's conserved quantities and environmental coupling, implemented via "fix" commands that integrate with time-stepping algorithms. The microcanonical (NVE) ensemble, corresponding to constant number of particles (N), volume (V), and energy (E), uses the fix nve command to perform velocity-Verlet integration on atom positions and velocities without external thermostats or barostats, preserving total energy for isolated systems. The canonical (NVT) ensemble maintains constant N, V, and temperature (T) using thermostats such as Nose-Hoover, implemented via fix nvt, which introduces an extended dynamical variable η to couple the system to a heat bath; the thermostat evolution follows the equation
\frac{d\eta}{dt} = \frac{T - T_\text{target}}{Q},
where T is the instantaneous temperature, T_target is the desired temperature, and Q is the thermostat mass parameter, typically set to allow damping over ~100 timesteps. Alternatively, the Langevin thermostat (fix langevin) adds stochastic forces and frictional drag to mimic viscous solvent effects, with damping parameter τ_damp controlling relaxation (e.g., 100 femtoseconds for water-like systems). The isothermal-isobaric (NPT) ensemble extends NVT by coupling to pressure (P) via barostats in fix npt, using Nose-Hoover chains to adjust box volume and shape anisotropically (aniso or tri keywords) while targeting a reference pressure, with damping times ~1000 timesteps for stability. Rigid body variants (e.g., fix rigid/nvt, fix rigid/npt) integrate these ensembles for constrained systems, accounting for removed degrees of freedom.[34][35][36][37][38]
Constraints in LAMMPS enforce physical restrictions on atomic motion, often to model molecular rigidity or confinement, applied through dedicated fix commands during dynamics or minimization. Bond, angle, and dihedral constraints are handled by fix shake or fix rattle, which iteratively solve for constraint forces to maintain specified lengths and angles (e.g., in water molecules), allowing larger timesteps (up to 2 fs); SHAKE adjusts positions via Lagrange multipliers, while RATTLE additionally constrains velocities to avoid drift, both converging within a tolerance like 10^{-4} over 20 iterations. Rigid body constraints via fix rigid or fix rigid/small treat atom clusters as indivisible units, computing center-of-mass forces and torques for collective motion, compatible with NVE, NVT, and NPT ensembles by adjusting thermodynamic degrees of freedom. Confinement is achieved with walls, such as planar (fix wall/reflect for elastic bounces or fix wall/lj93 for Lennard-Jones interactions) or spherical (fix wall/region with cylindrical/spherical regions), which bound particles and can move dynamically (e.g., oscillating or rotating) to simulate interfaces or pores.[39][38][40]
For non-equilibrium molecular dynamics (NEMD), LAMMPS enables transport property calculations by applying external drives, such as shear via fix deform to strain the box at controlled rates, combined with SLLOD integration in fix nvt/sllod to thermostat the peculiar velocity profile and compute viscosity from stress-strain relations. Thermal conductivity is assessed using fix thermal/conductivity for reverse NEMD, swapping velocities across slabs to impose heat flux. Monte Carlo methods supplement MD for reactive or compositional sampling; fix atom/swap performs atom-type exchanges with a reservoir at specified temperature, while fix bond/swap swaps bonds in polymers for equilibration, and fix gcmc enables grand canonical insertions/deletions of atoms or molecules, including reactions via template-based swaps in the REAX package. These integrate with MD timestepping for hybrid simulations, such as reactive force fields. As of the 22 July 2025 release, enhancements include fix hmc for Hamiltonian Monte Carlo and fix neighbor/swap for hybrid MD/kMC simulations.[41][42][43][44][45][46][47][24]
Time Integration and Algorithms
LAMMPS employs the velocity-Verlet algorithm as its primary time integration scheme for molecular dynamics simulations, implementing the velocity form of the Störmer-Verlet method to update atom positions and velocities while conserving energy in the microcanonical ensemble.[34] This integrator is invoked through fixes such asfix nve, which performs the updates during each timestep without additional thermodynamic manipulations.[34] The algorithm proceeds in sub-steps per timestep: an initial half-step velocity update, a full position update using the half-step velocity, force computation from the new positions, and a final half-step velocity update.[48]
The core updates in the velocity-Verlet scheme are given by the following equations, where \mathbf{r}(t) is the position at time t, \mathbf{v}(t) is the velocity, \mathbf{a}(t) is the acceleration derived from forces at time t, and \Delta t is the timestep size:
\mathbf{v}(t + \Delta t/2) = \mathbf{v}(t) + \frac{\Delta t}{2} \mathbf{a}(t)
\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \Delta t \cdot \mathbf{v}(t + \Delta t/2)
(Forces and accelerations are then recomputed at the new positions to obtain \mathbf{a}(t + \Delta t), followed by)
\mathbf{v}(t + \Delta t) = \mathbf{v}(t + \Delta t/2) + \frac{\Delta t}{2} \mathbf{a}(t + \Delta t)
These equations ensure time-reversible dynamics and second-order accuracy, making the method suitable for a wide range of atomic and molecular systems.[34] The full timestep cycle also involves clearing and recomputing forces from pairwise, bonded, and long-range interactions after the initial integration.[48] As of the 22 July 2025 release, enhancements include new fix styles for bosonic path integral molecular dynamics.[24]
For overdamped systems, LAMMPS supports Brownian dynamics through the fix langevin command paired with fix nve, which adds frictional drag and random forces to model implicit solvent effects, effectively neglecting inertial terms in Newton's equations.[37] The drag force is \mathbf{F}_f = -\gamma \mathbf{v}, where \gamma = m / \text{damp} and damp is a user-specified damping parameter, while the random force scales with \sqrt{k_B T m / (\text{damp} \Delta t)} to maintain the target temperature.[37] Rigid body integration is handled by fixes like fix rigid/nve, which apply specialized algorithms to treat clusters of atoms as rigid units, using methods such as those in Miller et al. (2002) for improved energy conservation in constant-volume simulations.[38]
Energy minimization in LAMMPS relaxes atomic configurations to local energy minima using several algorithms specified via the min_style command. The conjugate gradient (cg) method employs the Polak-Ribiere formulation, iteratively updating the search direction by combining the current force gradient with the previous direction, orthogonalizing when progress stalls.[49] For damped dynamics approaches, quickmin applies a damping parameter based on the projection of velocities onto forces, initializing velocities to zero for steepest descent-like relaxation.[49] The fire style extends this with a variable timestep and velocity rescaling to retain components parallel to forces, enhancing efficiency for complex energy landscapes, particularly in nudged elastic band calculations.[49]
To handle systems with disparate force timescales, LAMMPS implements the reversible reference-system propagator algorithm (rRESPA) via the run_style respa command, enabling multi-level hierarchical timestepping where inner loops use smaller timesteps for fast motions like bonds, and outer loops for slower interactions like non-bonded forces.[50] Loop factors define the ratio of timesteps across levels, with defaults placing bonded interactions at the innermost level and pairwise/kspace at the outermost; this is based on the formulation in Tuckerman et al. (1992).[50] Adaptive timestepping is available through fix dt/reset, which recomputes the timestep every N steps to cap the maximum atomic displacement at Xmax, using a conservative estimate from current velocities and forces, bounded by user-defined minimum and maximum values.[51]
The rerun mode, invoked by the rerun command, allows reprocessing of atom snapshots from dump files with updated force fields or parameters, effectively restarting dynamics without full simulation overhead by computing new forces on pre-stored trajectories.[52] This facilitates diagnostics or potential changes mid-analysis, reading snapshots sequentially and triggering zero-timestep runs for force evaluations.[52]
Atom Manipulation and Initialization
LAMMPS provides several mechanisms for initializing atomic configurations, allowing users to set up simulations with precise control over atom positions, types, and properties before dynamics commence. Initial setups can be achieved by reading pre-defined coordinates from input files or generating atoms programmatically within the simulation box. These methods ensure compatibility with various particle models, such as atomic or molecular systems.[53] The primary method for importing atom coordinates is theread_data command, which loads data from an ASCII or gzipped file in LAMMPS's standard data format. This file specifies the simulation box dimensions, atom types, and the "Atoms" section containing coordinates (x, y, z) along with properties like atom ID and type, depending on the chosen atom_style (e.g., atomic or charge). Coordinates outside the box are wrapped into it for periodic boundaries, enabling flexible initialization of complex structures. For formats like XYZ, external tools such as Atomsk are commonly used to convert files to LAMMPS-compatible data files before reading. The command supports multiple invocations with keywords like add append to incrementally build systems by appending atoms from additional files, adjusting IDs via offset, or shifting positions with shift.[54][55]
Atom creation without external files begins with the lattice command, which defines a periodic lattice structure as a foundation for subsequent atom placement. Supported styles include simple cubic (sc), body-centered cubic (bcc), face-centered cubic (fcc), hexagonal close-packed (hcp), and custom lattices, each with a specified scale factor representing the lattice constant in distance units or reduced density for Lennard-Jones units. For instance, lattice fcc 3.615 sets up an FCC lattice with a spacing of 3.615 Å, suitable for metals like copper. Keywords like orient allow rotation of lattice vectors, while basis positions multiple atom types within the unit cell, facilitating the initialization of multi-component crystals. This lattice is then used by other commands to populate the box.[56]
The create_atoms command leverages the defined lattice or other methods to insert atoms into specified regions of the simulation box. It supports styles such as box to fill the entire box with lattice atoms, region to target a geometric volume (e.g., a sphere or cylinder defined via the region command), random to place a specified number of atoms at random positions with overlap checks, or single for precise placement at given coordinates. For example, create_atoms 1 region disk adds type-1 atoms on the lattice within the "disk" region, with options like basis 1 0.0 0.0 0.0 to select specific unit cell positions. An overlap distance keyword prevents atom clashes during insertion, ensuring valid initial configurations. This command serves as an alternative to file-based input for generating uniform or patterned structures.[53]
To scale up small initial systems, the replicate command duplicates the current configuration by factors nx, ny, nz in each dimension, multiplying the atom count accordingly (e.g., 2x2x2 yields 8 times more atoms). It preserves properties like velocities and bonds, reassigning atom IDs and partitioning the enlarged box across processors. Keywords such as bbox optimize performance for large replications, making it ideal for extending periodic unit cells into larger simulation domains during setup.[57]
Post-creation manipulation refines the atomic arrangement using commands like displace_atoms, set, and delete_atoms. The displace_atoms command shifts a group of atoms via styles including move for uniform translation (e.g., displace_atoms all move 0.0 0.0 5.0), ramp for position-dependent displacement along a dimension, random for Gaussian perturbations with a seed, or rotate around a specified axis by an angle in degrees. Displacements are applied in lattice or box units and remap atoms into periodic boundaries, useful for introducing initial strains or randomizing lattice positions.[58]
The set command assigns or modifies per-atom attributes, such as type, charge, or mass, targeting selections by atom ID range (e.g., *), type (e.g., type 1-3), group, region, or molecule ID. For initialization, set type 2 changes all atoms to type 2, altering interaction potentials, while set charge * 1.0 assigns a uniform charge; variables enable spatially varying values (e.g., set charge * v_myvar). Fraction-based setting, like set type/fraction 2 0.5 12345, randomly assigns types to a subset of atoms, supporting probabilistic setups. This command requires communication across processors for consistency in parallel runs. As of the 12 June 2025 update, a new fix set command has been added for dynamic attribute modifications during simulation.[59][24]
Atom removal is handled by delete_atoms, which eliminates specified atoms to create voids or eliminate overlaps. Styles include group to delete all in a group-ID, region to target a volume (with mol yes for entire molecules), overlap to remove one atom from close pairs within a cutoff (e.g., 1.0 Å), random to delete a fraction or count of atoms with a seed, or variable for conditional removal. Keywords like compress yes reassigns IDs sequentially, and bond yes breaks associated bonds, ensuring clean configurations after deletion. For example, delete_atoms overlap 0.5 all all resolves initial overlaps in dense packings. As of the 10 September 2025 update, a "condense" option has been added to further compress atom IDs.[60][24]
For meshing purposes, LAMMPS includes the compute voronoi/atom command, which performs a Voronoi tessellation of atom positions using the Voro++ library to divide space into cells closer to each atom than others. It outputs per-atom arrays with cell volume (in distance cubed units), number of faces (indicating neighbors), and optionally surface area, restricted to a group or all atoms. Options like radius adjust for polydisperse particles via atom-style variables. This computation aids in analyzing or meshing atomic structures by providing geometric decompositions, though it requires the VORONOI package.[61]
Dynamic modifications during simulation allow ongoing atom addition via the fix deposit command, which inserts single atoms or molecules every M timesteps until N insertions, typically for deposition processes. Insertions occur in a specified region (e.g., above a surface), with near R ensuring minimum distance from existing atoms, and velocity ranges set via vxlo vxhi. For molecules, a template-ID defines the structure. Random seeds control placement attempts, with units in lattice or box coordinates; failed insertions trigger warnings but continue. An example is fix dep all deposit 1000 1 10 12345 [region](/page/Region) above near 1.0, adding 1000 type-1 atoms over time. This fix integrates with timestepping without affecting force computations.[62]
Outputs and Analysis
Data Output Formats
LAMMPS provides several mechanisms for exporting simulation data during runs, enabling users to capture progress, atomic configurations, and state information for later analysis or continuation. These outputs are controlled through input script commands, with flexibility in frequency, format, and content to accommodate diverse computational needs.[63] Log files serve as the primary record of simulation progress, outputting timestep-by-timestep summaries of key thermodynamic quantities such as total energy, kinetic energy, potential energy components (e.g., pairwise, bond, angle), temperature, pressure, and volume. By default, LAMMPS writes this information to both the screen and a file named "log.lammps" in a human-readable text format, with styles like "one" for a single line per timestep or "multi" for detailed multi-line breakdowns.[64][65] Thethermo command dictates the output frequency (e.g., every N timesteps), and customization is achieved via the thermo_style custom option, which allows inclusion of user-defined quantities derived from compute or fix commands, such as specific stress tensor components or custom variables.[64][66]
Dump files capture detailed snapshots of the atomic system at specified intervals, preserving per-atom data like positions, velocities, forces, types, charges, and image flags (for unwrapping coordinates across periodic boundaries). The dump command supports a variety of formats tailored to different visualization and analysis tools: text-based options include XYZ for simple coordinate dumps compatible with molecular viewers, ATOM for extended per-atom properties, custom for user-selected attributes, and extxyz (added as of April 2025) for extended XYZ format with additional metadata; binary formats like DCD and XTC offer compact storage for trajectories; HDF5 (via h5md style) enables hierarchical data organization with metadata; and VTK facilitates volume rendering of grid-based data.[67] Parallel I/O is integrated for scalability, using filename wildcards (e.g., "dump.%") to distribute output across processors, though not all formats (e.g., XYZ, VTK) support this natively.[67] Output frequency is adjustable via the every keyword or variables, allowing irregular intervals based on simulation conditions.[67] As of June 2025, LAMMPS also supports JSON format for structured input and output of simulation data.[68]
Restart files store the full simulation state in binary format, enabling seamless continuation of interrupted runs from exact timesteps. Generated by the restart command at multiples of a specified interval (e.g., every 1000 steps), these files include atom coordinates, velocities, forces, and simulation parameters like box dimensions and neighbor lists, but exclude certain dynamic elements (e.g., some fixes) that must be re-specified in the input script.[69][70] Parallel I/O is supported through wildcards in filenames, with options to control the number of files (e.g., one per processor group) for efficient writing on large-scale systems.[69]
Compression enhances efficiency for voluminous outputs, with LAMMPS supporting GZIP (via .gz suffix) for text-based dump and restart files, and Zstandard (Zstd, via .zstd) for high-ratio compression in modern formats.[67] These options reduce storage demands without altering the core data structure, and variable output frequency further optimizes by triggering writes only when needed, such as after equilibration phases. Custom thermodynamic dumps, often integrated with dump files, leverage compute and fix outputs for specialized summaries like per-atom energies or stress profiles.[64] Such formats are commonly used in post-processing workflows for trajectory analysis.[67]
Diagnostics and Post-Processing
LAMMPS provides a suite of diagnostic tools for computing quantities during simulations, enabling users to monitor and analyze system properties without interrupting the primary dynamics. These tools include computes for calculating instantaneous values such as per-atom or global metrics, and fixes for applying time-dependent modifications to the simulation. Additionally, introspection commands allow querying the current state of the simulation, while post-processing options facilitate analysis of saved data externally or through re-runs.[71][72][73] Computes in LAMMPS calculate and store global, per-atom, local, or per-grid quantities derived from the atom positions, velocities, and forces. Per-atom computes produce values for individual atoms, such as the stress tensor via thecompute stress/atom command, which computes the symmetric per-atom stress tensor with six components stored as a per-atom array. Coordination numbers can be determined using the compute coord/atom command, which counts neighboring atoms within a specified cutoff distance for each atom in a group, yielding a per-atom vector. Global computes provide system-wide properties; for example, the compute rdf command calculates the radial distribution function (RDF) and coordination number for a group of particles by binning pairwise distances. Another key global compute is compute msd, which estimates the mean-squared displacement (MSD) of atoms in a group, accounting for periodic boundaries:
\text{MSD} = \left\langle \left( \mathbf{r}(t) - \mathbf{r}(0) \right)^2 \right\rangle
where the average is taken over the atoms in the group and unwrapped trajectories to handle diffusion across boundaries.[74][75][76][77]
Fixes enable dynamic adjustments during the simulation, including time-dependent computations that modify or average quantities. The fix ave/time command performs time averaging of global, per-atom, or local quantities from specified inputs, such as computes or thermodynamic data, outputting averaged values at regular intervals to the screen or a file. For instance, it can average temperature or pressure over multiple timesteps to reduce noise in statistical data. The fix deform command applies straining to the simulation box, allowing controlled deformation such as uniaxial extension or shear, which is useful for studying mechanical responses under applied stress.[78][42]
Introspection in LAMMPS is supported through commands that query and report the internal state of the simulation. The thermo_style command customizes the output of thermodynamic information printed during the run, specifying quantities like energy, temperature, or pressure in formats such as "one" for single-line summaries or "multi" for detailed multi-line reports. The info command provides detailed diagnostics on the current LAMMPS process, including details on active computes, fixes, variables, and simulation parameters, aiding in debugging and validation of input scripts.[64][73]
Post-processing in LAMMPS extends analysis beyond the initial simulation run using external tools and built-in re-analysis capabilities. The Pizza.py toolkit, a Python 2-based collection developed at Sandia National Laboratories, offers utilities for visualization and data manipulation of LAMMPS outputs, including wrappers for integrating with molecular viewers like PyMOL for rendering atomic structures and VMD for real-time display of trajectories; a Python 3-compatible fork known as Pizza3 is available for modern environments.[79][80][81] The rerun command allows re-processing of dump files from a previous simulation without recomputing the dynamics, enabling the application of new computes or fixes to archived snapshots for efficient exploration of alternative analyses.[52]
Extensions and Integrations
Coupling with Other Software
LAMMPS supports integration with external software through multiple interfaces, enabling hybrid simulations that combine molecular dynamics with other computational methods. These couplings facilitate multiscale modeling, where LAMMPS handles atomic-scale interactions while interfacing with continuum or quantum mechanics codes. The primary mechanisms include library embedding, message-passing standards, scripting wrappers, and dynamic plugins, allowing flexible incorporation without extensive code modifications.[82] In library mode, LAMMPS can be compiled as a static or shared library and embedded into larger applications using C, C++, Fortran, or Python wrappers. This approach allows external codes to drive LAMMPS simulations by invoking its API functions, such as executing commands via strings or files and extracting simulation data like atom positions and forces. The C library API, defined inlibrary.h, uses opaque handles for LAMMPS instances and supports thread-safe operations when using appropriate MPI communicators. Fortran wrappers leverage ISO_C_BINDING for seamless integration, while C++ provides direct class access for advanced users. This mode is particularly useful for tightly coupled systems where LAMMPS serves as a computational engine within a host application.[83]
The Message Passing Driver Interface (MDI) provides a standardized protocol for loose coupling between LAMMPS and other codes, operating in a client-server model over MPI or sockets. LAMMPS can act as either a driver (client), sending requests for data like forces or energies, or an engine (server), responding to incoming queries. This enables hybrid simulations, such as quantum mechanics/molecular mechanics (QM/MM) setups, where LAMMPS computes classical forces on surrounding atoms while a quantum code handles a reactive region. Examples include direct interfaces with Quantum ESPRESSO for ab initio molecular dynamics and tools like PySCF or NWChem via Python wrappers included in the LAMMPS distribution. The MDI package in LAMMPS supports fixes like fix mdi/qm and fix mdi/qmmm to manage these interactions during timesteps.[84]
Python integration enhances LAMMPS's extensibility through the lammps module, which wraps the C library API using the ctypes module for dynamic loading. Users can instantiate a lammps object to run simulations from Python scripts, execute commands interactively, and access internal data via methods like extract_compute() for real-time analysis during runs. This supports scripting workflows, such as automating parameter sweeps or embedding LAMMPS within larger Python-based applications, and allows call-outs to external Python code mid-simulation using the python command. The module ensures compatibility with parallel execution via mpi4py, making it suitable for high-performance computing environments.[85]
Practical examples of these couplings include hybrid molecular dynamics-finite element method (MD-FEM) simulations, where atom positions from LAMMPS inform boundary conditions in a finite element solver, which in turn returns interpolated forces to LAMMPS atoms. Another application is multiscale modeling via MDI, such as coupling LAMMPS with quantum codes for concurrent atomistic-continuum simulations. Additionally, the VORONOI package demonstrates fix-based coupling by integrating the Voro++ library to compute Voronoi tessellations during LAMMPS runs.[82]
LAMMPS's plugin system allows users to extend functionality without recompiling the core executable by loading dynamic shared object (DSO) files at runtime. The plugin command supports loading, unloading, listing, and clearing plugins, which can add new styles for pairs, fixes, or commands. Plugins are loaded globally and persist across LAMMPS instances, with automatic discovery via the LAMMPS_PLUGIN_PATH environment variable. This mechanism is ideal for incorporating user-developed or third-party extensions, such as custom potentials, while maintaining binary compatibility with the PLUGIN package.[86]