VASP Ab Initio Molecular Dynamics (AIMD)

Run Born-Oppenheimer molecular dynamics with DFT forces at each step. Expensive but provides finite-temperature behavior, diffusion coefficients, and reaction dynamics.

When to Use

Thermal stability — check if a structure is stable at operating temperature
Diffusion — compute diffusion coefficients (e.g., Li-ion conductors)
Reaction dynamics — observe bond breaking/forming at finite temperature
Free energy sampling — metadynamics or thermodynamic integration
Amorphous structures — melt-quench to generate amorphous phases

Basic NVT Molecular Dynamics

from catgo.workflow import Workflow
from catgo.workflow.builtins import geo_opt, md

wf = Workflow("AIMD at 600K")
struct = wf.add_task("structure_input", structure=structure_json)

# Optional: relax first
opt = wf.add_task(geo_opt, structure=struct.output.structure,
                  system_name="relax")

# NVT MD at 600 K
run = wf.add_task(md, structure=opt.output.structure,
                  IBRION=0,      # Molecular dynamics
                  NSW=5000,      # Number of MD steps
                  POTIM=1.0,     # Time step in fs
                  TEBEG=600,     # Starting temperature (K)
                  TEEND=600,     # Ending temperature (K)
                  SMASS=0,       # Nose-Hoover thermostat (NVT)
                  ISIF=2,        # Fix cell shape/volume
                  system_name="AIMD_600K")

wf.submit()

MCP Workflow

catgo_workflow_engine(action="create", params={"name": "AIMD 600K"})

catgo_workflow_engine(action="add_task", params={
  "workflow_id": "wf_xxx",
  "task_type": "structure_input",
  "structure": "<json>"
})

catgo_workflow_engine(action="add_task", params={
  "workflow_id": "wf_xxx",
  "task_type": "md",
  "software": "vasp",
  "structure": "{{t_001.output.structure}}",
  "NSW": 5000,
  "POTIM": 1.0,
  "TEBEG": 600,
  "TEEND": 600,
  "SMASS": 0,
  "system_name": "AIMD_600K"
})

catgo_workflow_engine(action="submit", params={"workflow_id": "wf_xxx"})

Key Parameters

Parameter	Default	Purpose
IBRION	0	Molecular dynamics mode
NSW	1000	Number of MD steps
POTIM	1.0	Time step in femtoseconds
TEBEG	300	Initial temperature (K)
TEEND	300	Final temperature (K). Set equal to TEBEG for isothermal
SMASS	-1	Thermostat: -1=NVE, 0=Nose-Hoover NVT, >0=Nose mass
ISIF	2	Fix cell (NVT). Use ISIF=3 for NPT (rare in AIMD)
NBLOCK	1	Write trajectory every NBLOCK steps

Thermostat Selection (SMASS)

SMASS	Ensemble	Use case
-1	NVE (microcanonical)	Energy conservation test, short dynamics
0	NVT Nose-Hoover	Standard production MD at fixed T
1-3	NVT with Nose mass	Larger SMASS = slower T coupling (less perturbation)
-3	Langevin thermostat	Better T control for small systems

Recommendation: Use SMASS=0 (Nose-Hoover) for most production runs. Use SMASS=-1 (NVE) for energy conservation checks and very short equilibration diagnostics.

Temperature Ramp (Heating/Cooling)

To heat from 300 K to 1500 K (simulated annealing or melt-quench):

run = wf.add_task(md, structure=s,
                  NSW=10000,
                  POTIM=2.0,
                  TEBEG=300,     # Start at 300 K
                  TEEND=1500,    # Ramp to 1500 K
                  SMASS=0,
                  system_name="heating_ramp")

Time Step Selection (POTIM)

System	Recommended POTIM (fs)	Reason
Heavy elements (Pt, Au, Ru)	2.0	Heavy atoms, slow dynamics
Oxides (TiO2, RuO2)	1.0-1.5	O is light, moderate step needed
Light elements (H, Li)	0.5-1.0	Fast H vibrations need small step
Proton transfer	0.5	H requires fine time resolution

Rule of thumb: if the total energy drifts upward in NVE, reduce POTIM.

Performance Settings

AIMD is expensive. Optimize performance:

run = wf.add_task(md, structure=s,
                  NSW=5000,
                  POTIM=1.0,
                  TEBEG=600,
                  SMASS=0,
                  # Performance
                  ALGO="VeryFast",    # Fastest SCF convergence per step
                  NELM=60,            # Limit SCF steps (MD does not need tight SCF)
                  EDIFF=1e-4,         # Looser SCF for MD (still accurate forces)
                  LREAL="Auto",       # Real-space projection for speed
                  LWAVE=False,        # Do not write WAVECAR each step
                  NCORE=4,
                  system_name="AIMD")

EDIFF=1e-4 is acceptable for MD. Forces at 1e-4 SCF convergence are sufficiently accurate for MD trajectories. This saves 30-50% compute time vs 1e-5.

Supercell Size

AIMD requires large enough supercells to avoid finite-size effects:

Minimum: 64-100 atoms for bulk liquids/diffusion
Surfaces: use the slab supercell (typically 2x2 or 3x3 surface unit cell)
Small molecules on surface: existing slab supercell is usually fine

Build a supercell before MD:

catgo_structure(action="supercell", params={"scaling": [2, 2, 1]})

Output

The md task produces:

output.trajectory — atomic positions at each step (XDATCAR)
output.energy — energy vs time

Monitoring a Running MD

catgo_workflow_engine(action="status", params={"workflow_id": "wf_xxx"})
# Check NSW progress, temperature stability

catgo_analyze(action="convergence", params={"task_id": "t_md"})
# Energy vs step, temperature vs step

Troubleshooting

Problem	Fix
Temperature explodes	Reduce POTIM, check initial structure for overlapping atoms
Energy drift in NVE	Reduce POTIM, tighten EDIFF to 1e-5
SCF not converging at each step	Use ALGO=VeryFast, increase NELM
Too slow	Reduce ENCUT (400 eV ok for MD), use LREAL=Auto, fewer k-points
Atoms evaporating from slab	Add more vacuum, or constrain bottom layers

Typical Simulation Lengths

Purpose	NSW	POTIM	Total time
Quick stability check	1000	1.0	1 ps
Equilibration	5000	1.0	5 ps
Production (diffusion)	20000-50000	1.0-2.0	20-100 ps
Melt-quench	10000+	2.0	20+ ps

vasp-md

SKILL.md.claude/skills/vasp-md/SKILL.md