How to Validate Molecular Dynamics Simulation Results Before Publishing
Running an MD simulation is the beginning, not the end. Before any result from a trajectory can be reported, it must pass a series of validation checks — not optional best practices, but the minimum standard that reviewers at major journals expect to see documented in your methods section.
Why validation is not optional
MD simulation produces output regardless of whether the simulation was set up correctly. A system with a non-zero charge, mismatched force fields, or insufficient equilibration runs to completion, generates a trajectory file, and returns numbers that look like results. The only way to distinguish a valid simulation from an invalid one is to run the checks described in this article.
This matters practically because reviewers at journals like the Journal of Chemical Theory and Computation, Journal of Chemical Information and Modeling, Journal of the American Chemical Society, and PLOS Computational Biology increasingly request validation data explicitly. Papers without it are sent back for revision. Papers that falsely imply validation and are then scrutinized post-publication create reputational damage that is difficult to recover from.
The five checks below are organized in the order you should run them — from basic system integrity through convergence. Work through them sequentially after every simulation before writing up any results.
Check 1 — Energy conservation and stability
Extract and plot temperature, pressure, potential energy, and total energy over the entire simulation (including equilibration phases) using gmx energy:
# Extract key thermodynamic properties
echo "Temperature Pressure Potential Total-Energy Density" | \
gmx energy -f md.edr -o thermo.xvg
# Plot with any tool — xmgrace, python/matplotlib, or gnuplotFor the production phase specifically, check for energy drift using the gmx check utility — a drift rate above 0.1 kJ/mol/atom/ns warrants investigation:
gmx check -e md.edr 2>&1 | grep -i "drift\|conservation"Check 2 — Equilibration confirmation with RMSD
Calculate backbone RMSD over the full trajectory — equilibration phases included. The production phase begins only when RMSD has stopped rising and plateaued into a stable range. Any data collected before plateau is not from an equilibrated system and must be excluded from analysis.
# Calculate RMSD over entire run including equilibration
gmx rms \
-s md.tpr \
-f md_centered.xtc \
-o rmsd_full.xvg \
-tu ns
# Select Backbone for both groupsCheck 3 — Structural integrity
RMSD tells you how much the structure moved. Structural integrity checks tell you whether it remained properly folded. Three complementary analyses cover this:
Radius of gyration — overall compactness
gmx gyrate -s md.tpr -f md_centered.xtc -o gyrate.xvg
# Select Protein
# A stable Rg means the protein maintained its overall foldSecondary structure assignment over time
GROMACS can track secondary structure (helix, sheet, coil) at every frame using the DSSP algorithm:
gmx do_dssp \
-s md.tpr \
-f md_centered.xtc \
-o ss.xpm \
-sc ss_count.xvg
# Select Protein
# Convert XPM to image for visual inspection
gmx xpm2ps -f ss.xpm -o ss.epsA protein that maintains its fold shows consistent secondary structure content throughout the simulation — helices and sheets at the same residue positions as in the starting structure, with only minor fluctuations at loop regions. Sudden loss of secondary structure elements or persistent coil in regions that should be helical is a flag for unfolding or force field problems.
Compare RMSF to experimental B-factors
If an experimental crystal structure is available, compare the per-residue flexibility from your simulation (RMSF) to the B-factors from crystallography. Strong correlation — flexible loops in the simulation matching high B-factor regions in the crystal — validates that your simulation is capturing real protein dynamics rather than force field artifacts:
# Generate RMSF with B-factor equivalent output
gmx rmsf \
-s md.tpr \
-f md_centered.xtc \
-oq bfactors.pdb \
-res
# bfactors.pdb has RMSF-converted B-factors in the B-column
# Compare visually in PyMOL: spectrum b, blue_red, bfactorsCheck 4 — Convergence across the trajectory
Equilibration checks confirm the system reached a stable state. Convergence checks confirm the system stayed there and that your observations are representative of equilibrium rather than one particular trajectory segment.
The standard convergence test is the block averaging or first-half vs second-half comparison: split the production trajectory at the midpoint and compare key metrics between the two halves.
# Split trajectory at midpoint for convergence check
# Assuming 100 ns production: first half = 0-50 ns, second = 50-100 ns
gmx rmsf \
-s md.tpr \
-f md_centered.xtc \
-o rmsf_first_half.xvg \
-b 0 -e 50000 \
-res
gmx rmsf \
-s md.tpr \
-f md_centered.xtc \
-o rmsf_second_half.xvg \
-b 50000 -e 100000 \
-res
# Compare the two profiles — they should be similar
# Systematic differences indicate the simulation hasn't convergedCheck 5 — Independent replica agreement
A single MD trajectory is a single sample from a stochastic process. Key conclusions — especially quantitative ones like binding free energies, conformational populations, or H-bond occupancies — should be reproduced across at least two independent simulations started from different initial velocities. This is the strongest validation you can provide.
Run a second replica using the same preparation and simulation settings, but generate different initial velocities by using a different random seed in grompp:
# Generate replica with different random seed
gmx grompp \
-f md.mdp \
-c npt.gro \
-t npt.cpt \
-p topol.top \
-o md_rep2.tpr \
-seed 87654321 # different seed = different initial velocities
gmx mdrun -v -deffnm md_rep2Compare RMSD, RMSF, key distances, and for protein-ligand systems, ligand RMSD and H-bond occupancy between replicas. Conclusions that appear in both replicas are robust; those that appear in only one are either artifacts or require much longer sampling to converge.
What to include in your methods section
A reproducible, publication-ready methods section for an MD study should contain every parameter and decision needed for another researcher to replicate your work. Here is an annotated example of what reviewers at computational chemistry journals expect:
Every highlighted element in that paragraph is a specific detail that reviewers look for. Missing any of them is likely to generate a revision request asking you to supply it.
- Temperature, pressure, potential energy plotted — all stable throughout production
- RMSD plotted over full trajectory — production phase begins after clear plateau
- Analysis window explicitly defined and limited to post-equilibration frames
- Rg stable — protein maintained its fold
- Secondary structure content consistent throughout production
- First-half vs second-half RMSF comparison passes — similar profiles
- At least two independent replicas run — key conclusions reproduced in both
- Methods section includes: software version, force fields, box type and size, water model, salt concentration, equilibration protocol, timestep, constraints, electrostatics cutoff, production length, analysis window
Validation in five checks
Energy conservation confirms the simulation ran physically correctly. RMSD plateau confirms the system equilibrated before data collection. Structural integrity checks confirm the protein maintained its fold. Convergence checks confirm observations are reproducible across the trajectory. Independent replicas confirm conclusions are not artifacts of a single stochastic path. Run all five before reporting any result. Document all five in your methods section. Reviewers have seen every shortcut — the ones who read your methods are specifically checking for each of these.