Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Practical Molecular Dynamics with LAMMPS

Objective: Understand the practical implementation of molecular dynamics simulations, from theory to executable code, using LAMMPS as the implementation framework.

1. Introduction: From Theory to Practice

In previous lectures, we’ve covered the theoretical foundations of molecular dynamics:

Now we must translate this theory into actual simulations. This requires understanding:

  1. How MD codes are structured

  2. What decisions you must make as the researcher

  3. How to interpret and validate results

2. The MD Software Landscape

2.1 Major MD Packages

SoftwareStrengthsTypical Applications
LAMMPSMaterials, generality, parallelizationMetals, polymers, granular materials, simple fluids
GROMACSBiomolecules, speed, analysis toolsProteins, lipids, drug discovery
NAMDBiomolecules, GPU accelerationLarge biomolecular systems, membrane proteins
AMBERForce fields for biomoleculesDNA, RNA, proteins
VASP/CP2KAb initio MD (quantum mechanics)Reactive chemistry, electronic structure

Why LAMMPS for this course?

2.2 Common Architecture

Despite differences, most MD codes share a common structure:

┌─────────────────────────────────────┐
│   1. INITIALIZATION                 │
│   - Set units, boundary conditions  │
│   - Define atom types               │
└──────────────┬──────────────────────┘
               ↓
┌─────────────────────────────────────┐
│   2. SYSTEM BUILDING                │
│   - Create simulation box           │
│   - Add atoms (from file or lattice)│
│   - Set initial velocities          │
└──────────────┬──────────────────────┘
               ↓
┌─────────────────────────────────────┐
│   3. FORCE FIELD                    │
│   - Define interactions (LJ, bonds) │
│   - Set cutoffs and parameters      │
└──────────────┬──────────────────────┘
               ↓
┌─────────────────────────────────────┐
│   4. SIMULATION PROTOCOL            │
│   - Set timestep                    │
│   - Choose ensemble (NVE/NVT/NPT)   │
│   - Define output frequency         │
└──────────────┬──────────────────────┘
               ↓
┌─────────────────────────────────────┐
│   5. EXECUTION                      │
│   - Energy minimization (optional)  │
│   - Equilibration run               │
│   - Production run                  │
└──────────────┬──────────────────────┘
               ↓
┌─────────────────────────────────────┐
│   6. ANALYSIS                       │
│   - Extract thermodynamic averages  │
│   - Calculate structural properties │
│   - Visualize trajectories          │
└─────────────────────────────────────┘

3. The LAMMPS Input Script Philosophy

3.1 Sequential Execution

Unlike compiled programs where you define functions and then call them, LAMMPS executes commands line-by-line, top to bottom. Think of it as a recipe:

Step 1: Preheat oven to 350°F    →  units lj
Step 2: Mix flour and sugar      →  create_atoms 1 box
Step 3: Add eggs                 →  velocity all create 1.0 12345
Step 4: Bake for 30 minutes      →  run 50000

This means:

3.2 The Command Structure

Every LAMMPS command has a specific syntax:

command_name  argument1  argument2  keyword1 value1  keyword2 value2

Example:

fix  1  all  nvt  temp 300.0 300.0  100.0
│    │   │    │    │         │       │
│    │   │    │    └─────────┴───────┴─ Required arguments
│    │   │    └───────────────────────── Fix style (Nosé-Hoover)
│    │   └────────────────────────────── Atom group to apply fix
│    └────────────────────────────────── Fix ID (user-defined name)
└─────────────────────────────────────── Command name

3.3 Comments and Formatting

# This is a comment (ignored by LAMMPS)

units lj                    # Inline comments are allowed

# Blank lines are ignored - use them to organize sections

# Long commands can be split with &
pair_coeff  1 1  1.0 1.0 &
            2.5              # Cutoff on next line

4. Critical Decision Points

4.1 Unit Systems

LAMMPS supports multiple unit systems. This is the #1 source of errors.

Unit SystemLengthTimeEnergyMassTemperature
ljστεmkBT/ϵk_B T/\epsilon
realÅfskcal/molamuK
metalÅpseVamuK
simsJkgK

Golden Rule: All parameters in your input file must be consistent with the chosen unit system.

Example Error:

units real                           # Using real units (kcal/mol)
pair_coeff 1 1  1.0 1.0 2.5          # But these are LJ reduced units!
                                     # ✗ WRONG - Will give nonsense results

Correct:

units real
pair_coeff 1 1  0.238 3.405 10.0     # Argon in real units (kcal/mol, Å)

4.2 Reduced vs. Real Units

For simple systems (noble gases, generic polymers), reduced LJ units are preferred:

For specific materials, use real units:

Conversion Example (Argon):

In reduced units:

units lj
pair_coeff 1 1  1.0 1.0 3.0          # ε=1, σ=1, cutoff=3σ
velocity all create 1.0 12345        # T* = 1.0

In real units:

units real
pair_coeff 1 1  0.238 3.405 10.215   # ε in kcal/mol, σ in Å
velocity all create 119.8 12345      # T = 119.8 K (same as T*=1.0)

Where:

5. Ensemble Selection

5.1 NVE: The Fundamental Ensemble

Microcanonical Ensemble — Constant N, V, E

timestep 0.005
run 100000                           # No thermostat/barostat = NVE

Uses:

What to check:

5.2 NVT: Constant Temperature

Canonical Ensemble — Constant N, V, T

fix 1 all nvt temp 300.0 300.0 100.0

Parameters:

When to use:

Alternatives:

fix 1 all langevin 300.0 300.0 100.0 12345    # Langevin dynamics
fix 1 all temp/berendsen 300.0 300.0 100.0    # Berendsen (equilibration only!)

5.3 NPT: Constant Pressure

Isothermal-Isobaric Ensemble — Constant N, P, T

fix 1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0

Parameters:

Pressure control modes:

When to use:

Critical:

6. Timestep Selection

6.1 The Stability Criterion

The timestep Δt\Delta t must be much smaller than the fastest motion in the system.

Physical basis: Atomic vibrations have a period TvibT_{vib}. To accurately resolve this motion:

Δt≪Tvib\Delta t \ll T_{vib}
(1)

Typically: Δt≈Tvib/10\Delta t \approx T_{vib}/10

6.2 Timestep Guidelines

System TypeFastest MotionTypical Δt\Delta t
Argon (no bonds)LJ vibrations5 fs (real), 0.005 (LJ)
Water (rigid bonds)Constrained H-O2 fs
Water (flexible bonds)O-H stretch0.5 fs
Proteins (SHAKE)Constrained bonds2 fs
Proteins (flexible)C-H stretch1 fs
Coarse-grainedSoft potentials10-50 fs

In LAMMPS:

# Reduced LJ units
timestep 0.005                       # 0.005 τ

# Real units (femtoseconds)
timestep 1.0                         # 1 fs

6.3 Testing Your Timestep

Method: Run short NVE simulations with different timesteps and plot total energy.

# Test run
timestep 0.001
fix 1 all nve
run 10000
write_data test_0.001.data

# Repeat with 0.005, 0.01, 0.02...

Good timestep: Total energy fluctuates around a constant value (±0.01%) Too large: Total energy drifts upward (exponentially growing errors)

7. Output and Data Collection

7.1 Thermodynamic Output

thermo_style custom step temp press pe ke etotal vol density
thermo 100                           # Print every 100 steps

Column selection:

Writing to file:

log thermo.out                       # Redirect output to file

7.2 Trajectory Output (Snapshots)

dump 1 all custom 1000 traj.lammpstrj id type x y z vx vy vz
dump_modify 1 sort id                # Sort by atom ID for consistency

Frequency considerations:

File formats:

7.3 Restart Files

restart 50000 restart.*.data         # Save restart every 50k steps
write_restart final.restart          # Save at end

Use cases:

Restarting:

read_restart final.restart
# Continue with new commands
run 100000

8. The Multi-Stage Simulation Protocol

Real research simulations are rarely a single run command. The standard workflow has multiple stages:

8.1 Stage 1: Energy Minimization

Purpose: Remove steric clashes (overlapping atoms)

minimize 1.0e-4 1.0e-6 1000 10000
#        etol   ftol   maxiter maxeval

When to use:

Output to check:

Minimization converged in 547 steps
Final energy: -8342.3 kcal/mol

If minimization fails (doesn’t converge):

8.2 Stage 2: NVT Equilibration

Purpose: Heat system to target temperature

# Start from 0 K (minimized structure has no velocities)
velocity all create 0.0 12345

# Ramp temperature from 0 to 300 K over 50 ps
fix 1 all nvt temp 0.0 300.0 100.0
run 10000                            # 50 ps at dt=0.005 ps
unfix 1                              # Remove this thermostat

Duration: 10,000-100,000 steps depending on system size

8.3 Stage 3: NPT Equilibration (if needed)

Purpose: Equilibrate density at target pressure

fix 1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run 100000                           # 500 ps - pressure is slow!
unfix 1

What to monitor:

8.4 Stage 4: Production Run

Purpose: Collect data for analysis

# Switch to final ensemble (e.g., NVT for canonical averages)
fix 1 all nvt temp 300.0 300.0 100.0

# Reset counters and start fresh output
reset_timestep 0
log production.out

# Long run for statistics
run 1000000                          # 5 ns

write_data final_production.data

Duration: System-dependent

9. Practical Considerations for HPC

9.1 Parallel Efficiency

LAMMPS uses spatial decomposition — the simulation box is divided among processors.

Optimal scaling:

Beyond this: Communication overhead dominates (diminishing returns)

Test scaling:

# Run with different core counts
mpirun -np 4 lmp -in input.in        # Record walltime
mpirun -np 8 lmp -in input.in
mpirun -np 16 lmp -in input.in
# Ideal: 2x cores = 2x speedup (rarely achieved)

9.2 I/O Optimization

Disk writes are slow — they can dominate runtime for large systems.

Strategies:

Example:

# During equilibration - minimal output
thermo 10000
dump 1 all custom 50000 eq.lammpstrj id type x y z

# During production - detailed output
thermo 1000
dump 2 all dcd 1000 prod.dcd

9.3 Memory Management

LAMMPS memory scales with:

Typical usage: ~100-200 bytes/atom

For 1 million atoms: ~200 MB RAM (plus overhead)

If you run out of memory:

10. Common Pitfalls and Debugging

10.1 Energy Explosion

Symptom: Total energy increases exponentially, simulation crashes

Causes:

  1. Timestep too large → Reduce by factor of 2

  2. Bad initial configuration (overlapping atoms) → Run energy minimization

  3. Wrong units (mixing real and LJ) → Check all parameters

  4. Force field mismatch → Verify pair_coeff values

Diagnostic:

# Quick test with very small timestep
timestep 0.0001
run 100

If this works, gradually increase timestep.

10.2 System Doesn’t Equilibrate

Symptom: Temperature/pressure/density still drifting after 100,000 steps

Causes:

  1. Equilibration time too short → Run 10× longer

  2. Multiple metastable states (glass) → Need enhanced sampling

  3. Box too small → Finite-size effects

  4. Wrong ensemble (using NVE when you meant NVT)

Diagnostic: Plot running average of property — if slope ≠ 0, not equilibrated.

10.3 Nonsensical Results

Symptom: Calculated properties don’t match literature or expectations

Checklist:

11. Validation and Verification

11.1 Internal Consistency Checks

Before trusting any result:

  1. Energy conservation (NVE test):

    fix 1 all nve
run 100000
# Check: Total energy drift < 0.01% per ns

  2. Temperature distribution (NVT test):

    • Calculate kinetic energy distribution

    • Should follow chi-squared with NfN_f degrees of freedom

    • Tutorial 7 shows how to check this

  3. Pressure/volume stability (NPT test):

    • Running average of volume should be flat

    • Pressure should fluctuate around target (±10-20% is normal)

11.2 Comparison with Known Results

For learning/testing:

Example: Argon liquid at 94.4 K

Literature values (from NIST):

Your simulation should match within 1-2%.

12. From Simulation to Publication

12.1 What to Report

Minimum information for reproducibility:

  1. System details:

    • Number of atoms, composition

    • Box size (initial and final if NPT)

    • Initial configuration (lattice, random, from file)

  2. Force field:

    • Potential type (LJ, AMBER, etc.)

    • All parameters (ε, σ, charges, bonds)

    • Cutoffs and long-range corrections

  3. Simulation protocol:

    • Timestep

    • Ensemble (NVE/NVT/NPT)

    • Thermostat/barostat type and parameters

    • Equilibration duration

    • Production duration

  4. Analysis method:

    • How equilibration was detected

    • How averages were calculated

    • Uncertainty estimation method (block averaging)

Example (Methods section):

We performed molecular dynamics simulations of 2000 Argon atoms using the Lennard-Jones potential with ε = 0.238 kcal/mol and σ = 3.405 Å, with a cutoff of 10.0 Å. Simulations were performed using LAMMPS version 20210310 on the Purdue Anvil cluster. The system was equilibrated in the NPT ensemble at T = 94.4 K and P = 1 atm for 100 ps, followed by a 1 ns production run in the NVT ensemble. The Nosé-Hoover thermostat with damping parameter 100 fs was used for temperature control. A timestep of 1 fs was used throughout. Average density was calculated by discarding the first 200 ps of the production run and block-averaging the remaining 800 ps with 10 blocks, yielding ρ = 1.372 ± 0.008 g/cm³ (95% CI).

12.2 Module Loading for Reproducibility

Always document the exact modules used:

# Load required modules (Anvil)
module load gcc/11.2.0 openmpi/4.0.6
module load lammps/20210310

This ensures others can reproduce your results with the same software environment.

12.3 Data Archiving

Make your simulations reproducible:

13. Summary: The Simulation Checklist

Before running any production simulation:

Planning:

Setup:

Execution:

Analysis:

Reporting:

14. Further Reading

14.1 Essential References

Textbooks:

LAMMPS Documentation:

Papers on Best Practices:

14.2 Online Resources

15. Looking Ahead

This lecture covered the core workflow of running MD simulations. Advanced topics (covered in future lectures or independent study):

The skills you’ve learned — writing input files, running simulations on HPC, analyzing output — form the foundation for all of these advanced techniques.

Lecture Notes
Introduction to Molecular Dynamics: Concepts, Algorithms, and Formalism
Lecture Notes
Classical Interatomic Potentials: From Pair Interactions to Molecular Force Fields