GROMACS Tutorial for Beginners: Run Your First MD Simulation From Scratch (2026)
By the end of this tutorial you will have run a complete molecular dynamics simulation — protein prepared, equilibrated, production MD completed, and basic analysis done. Every command is exact and explained. No prior MD experience required.
What you need before starting
This tutorial continues from where the protein preparation tutorial ends. Before starting, you should have:
- GROMACS installed and working —
gmx --versionreturns a version number - A prepared, energy-minimized system:
em.gro,topol.top, andposre.itp - If you haven’t done the preparation yet, follow the protein preparation tutorial first and come back here
Force field: CHARMM36 with SPC/E water.
Simulation length: 1 ns production MD — short enough to finish in under an hour on a modern GPU, or a few hours on CPU.
Working directory structure
All commands in this tutorial run from ~/md/lysozyme/. At the start of this tutorial your directory should contain:
ls ~/md/lysozyme/
# You should see: em.gro em.edr em.log topol.top posre.itp ions.gro ...The most important files going forward are em.gro (your minimized coordinates), topol.top (the complete topology), and posre.itp (position restraints used during equilibration).
Steps 1–6: System preparation recap
If you followed the preparation tutorial, your system is already minimized and ready. Here is a one-command summary of where you are:
✓ Topology generated with
pdb2gmx — CHARMM36, SPC/E water✓ Box defined: dodecahedron, 1.0 nm minimum distance
✓ System solvated: ~7,000 water molecules
✓ Ions added: system neutral at 0.15 M NaCl
✓ Energy minimization converged: Fmax < 1000 kJ/mol/nm
Step 7 — NVT equilibration
NVT equilibration heats the system to its target temperature (300 K) while holding the volume fixed. The protein is restrained with position restraints during this phase — the water and ions equilibrate thermally around the protein without the protein itself moving dramatically. This typically runs for 100 ps.
Create the NVT parameter file:
; Run control
integrator = md ; Leap-frog integrator
nsteps = 50000 ; 50,000 steps × 2 fs = 100 ps
dt = 0.002 ; 2 femtosecond timestep
; Output
nstxout = 500 ; Save coordinates every 1 ps
nstvout = 500
nstenergy = 500 ; Save energies every 1 ps
nstlog = 500
; Bond constraints
constraint_algorithm = lincs ; LINCS — faster than SHAKE
constraints = h-bonds ; Constrain X-H bonds — enables 2 fs timestep
lincs_iter = 1
lincs_order = 4
; Non-bonded interactions
cutoff-scheme = Verlet
ns_type = grid
nstlist = 10
rcoulomb = 1.2 ; Electrostatics cutoff (nm)
rvdw = 1.2 ; Van der Waals cutoff (nm)
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4
fourierspacing = 0.16
; Temperature coupling — V-rescale thermostat
tcoupl = V-rescale ; Correct canonical ensemble
tc-grps = Protein Non-Protein ; Couple protein and solvent separately
tau_t = 0.1 0.1 ; Temperature coupling time constant (ps)
ref_t = 300 300 ; Target temperature 300 K for both groups
; Pressure — not controlled in NVT
pcoupl = no
; Periodic boundary conditions
pbc = xyz
; Position restraints — defined in posre.itp
define = -DPOSRES ; Apply position restraints on protein heavy atoms
# Pre-process
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
# Run NVT equilibration
gmx mdrun -v -deffnm nvtThis runs for 100 ps and takes a few minutes on CPU, under a minute on GPU. When complete, check that temperature stabilized at 300 K:
echo "Temperature" | gmx energy -f nvt.edr -o temperature.xvg
# Plot temperature.xvg — should plateau at ~300 Kdefine = -DPOSRES line activates the position restraints defined in posre.itp. These apply a harmonic force to all protein heavy atoms, keeping them close to their starting positions while the water and ions equilibrate freely around them. Without this, the protein can drift significantly during the early violent phase of heating.
Step 8 — NPT equilibration
NPT equilibration adds pressure control to allow the system density to reach its equilibrium value at 1 bar. The position restraints remain active. This phase also runs for 100 ps. The NPT MDP file is almost identical to NVT — add the pressure coupling block and remove the pcoupl = no line:
; Pressure coupling — Parrinello-Rahman barostat
pcoupl = Parrinello-Rahman ; Correct NPT ensemble
pcoupltype = isotropic ; Isotropic box scaling
tau_p = 2.0 ; Pressure coupling time constant (ps)
ref_p = 1.0 ; Target pressure 1.0 bar
compressibility = 4.5e-5 ; Compressibility of water (bar⁻¹)
refcoord_scaling= com ; Scale center of mass with box
; Keep position restraints active
define = -DPOSRES
# Pre-process (use nvt.gro as starting coordinates)
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
# Run NPT equilibration
gmx mdrun -v -deffnm nptThe -t nvt.cpt flag passes the checkpoint file from the NVT run, which contains the final velocities — you’re continuing the simulation with the same kinetic energy rather than randomly reinitializing it. Always pass checkpoint files between consecutive simulation phases.
After NPT completes, verify both temperature and pressure stabilized:
echo "Pressure" | gmx energy -f npt.edr -o pressure.xvg
echo "Density" | gmx energy -f npt.edr -o density.xvg
# Pressure: fluctuates around 1 bar (large fluctuations are normal)
# Density: should stabilize near 1000 kg/m³Step 9 — Production MD
Production MD removes all position restraints and runs the actual simulation. The key differences from equilibration: no -DPOSRES define, the thermostat switches to Nosé-Hoover for the correct canonical ensemble, and the run is longer. For this tutorial we run 1 nanosecond (500,000 steps at 2 fs).
integrator = md
nsteps = 500000 ; 500,000 × 2 fs = 1 ns
dt = 0.002
; Output — save every 10 ps
nstxout-compressed = 5000 ; Compressed trajectory (xtc) — saves disk space
nstenergy = 5000
nstlog = 5000
; Constraints — same as equilibration
constraint_algorithm = lincs
constraints = h-bonds
; Non-bonded — same as equilibration
cutoff-scheme = Verlet
nstlist = 10
rcoulomb = 1.2
rvdw = 1.2
coulombtype = PME
; Temperature — Nosé-Hoover for correct ensemble in production
tcoupl = Nose-Hoover
tc-grps = Protein Non-Protein
tau_t = 0.1 0.1
ref_t = 300 300
; Pressure — Parrinello-Rahman
pcoupl = Parrinello-Rahman
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 4.5e-5
; No position restraints in production
pbc = xyz
DispCorr = EnerPres ; Dispersion correction for CHARMM36
# Pre-process using NPT final coordinates and checkpoint
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
# Run production — GPU acceleration used automatically if available
gmx mdrun -v -deffnm mdExpected runtime: approximately 30–90 minutes on a GPU, 4–8 hours on CPU for a 1 ns simulation of lysozyme. GROMACS prints estimated completion time after the first few steps. For longer simulations, use -ntmpi 1 -ntomp 8 to specify thread counts on multi-core hardware.
grompp is fast), copy the md.tpr to the cluster, and submit gmx mdrun -v -deffnm md -ntmpi 1 -ntomp 8 -gpu_id 0 in your job script. The .tpr file is portable across machines and GROMACS versions.
Step 10 — Basic analysis
First: center and wrap the trajectory
Because of periodic boundary conditions, the protein may appear fragmented or drift across the box boundary in the raw trajectory. Always process it first:
gmx trjconv \
-s md.tpr \
-f md.xtc \
-o md_centered.xtc \
-center \
-pbc mol \
-ur compactWhen prompted, select Protein for centering and System for output. This produces a clean trajectory where the protein stays centered in the box throughout.
Calculate RMSD
RMSD measures how much the protein structure has deviated from the starting structure over time. A plateau indicates the protein has reached a stable conformation.
gmx rms \
-s md.tpr \
-f md_centered.xtc \
-o rmsd.xvg \
-tu nsSelect Backbone for both the least-squares fit group and the RMSD group. The -tu ns flag converts the x-axis to nanoseconds. Plot the resulting XVG file — for a well-behaved 1 ns lysozyme simulation, RMSD typically plateaus around 0.1–0.2 nm after 0.2–0.3 ns.
0.000 0.0000
0.010 0.0523
0.100 0.1187
0.500 0.1641
1.000 0.1728
RMSD plateaus around 0.17 nm — the protein backbone is stable throughout the 1 ns simulation.
Calculate RMSF (per-residue flexibility)
gmx rmsf \
-s md.tpr \
-f md_centered.xtc \
-o rmsf.xvg \
-resSelect Backbone. The -res flag averages over residues rather than atoms, giving you one flexibility value per residue. High RMSF residues are flexible loops; low RMSF residues are rigid helices and sheets. Plotting this against residue number gives you the protein’s flexibility fingerprint.
Calculate radius of gyration
gmx gyrate -s md.tpr -f md_centered.xtc -o gyrate.xvgSelect Protein. Radius of gyration measures overall compactness — a stable protein shows a near-constant Rg throughout the simulation. Decreasing Rg indicates compaction; increasing Rg indicates unfolding or domain separation.
All analysis tools at a glance
| GROMACS tool | What it calculates | Input group to select |
|---|---|---|
| gmx rms | RMSD vs reference structure — overall structural drift over time | Backbone (fit + RMSD) |
| gmx rmsf | Per-residue flexibility — which regions are rigid vs flexible | Backbone (+ -res flag) |
| gmx gyrate | Radius of gyration — overall protein compactness | Protein |
| gmx hbond | Hydrogen bond count over time — intra- and inter-molecular | Protein / Protein-Ligand |
| gmx energy | Thermodynamic properties — temperature, pressure, potential energy | Select from numbered list |
| gmx sasa | Solvent-accessible surface area — protein exposure to solvent | Protein |
What to do next
You’ve completed a full MD simulation pipeline — from raw PDB to trajectory to analysis. Here is how to build on what you’ve learned:
- Run longer: Increase
nstepsinmd.mdpto 50,000,000 (100 ns) for production-quality results - Add a ligand: The protein-ligand MD tutorial on this site extends this workflow to simulate drug-protein complexes
- Improve analysis: The trajectory analysis tutorial covers RMSD clustering, hydrogen bond occupancy, and MM-GBSA calculation
- Visualize the trajectory: Load
md_centered.xtcandmd.tprinto VMD or PyMOL to watch the protein move
- NVT equilibration completed — temperature stable at 300 K
- NPT equilibration completed — density stable near 1000 kg/m³
- Production MD completed — 1 ns trajectory saved as md.xtc
- Trajectory centered and wrapped with gmx trjconv
- RMSD calculated — backbone RMSD plateaus during production
- RMSF calculated — per-residue flexibility profile generated
You’ve run your first MD simulation
The workflow you’ve completed — preparation, minimization, NVT, NPT, production, analysis — is the same sequence used in every all-atom protein MD publication. The protein changes; the pipeline doesn’t. Swap 1AKI for any protein PDB ID, repeat the preparation tutorial, and this exact set of MDP files will run your simulation. Every parameter choice here is a reasonable default — the next step is understanding what each parameter does so you can adjust it when your research question demands it.