10 Common Molecular Dynamics Simulation Mistakes (And How to Avoid Every One)

Periodic boundary conditions require the simulation box to be large enough that the protein never comes within interaction distance of its own periodic image. When the minimum distance between the protein and the box edge is smaller than the non-bonded cutoff (typically 1.0–1.2 nm), the protein effectively interacts with itself — an unphysical artifact that corrupts electrostatics, introduces false forces, and produces meaningless dynamics.

This is especially problematic for elongated proteins, flexible proteins that sample extended conformations during the simulation, and any system that undergoes significant conformational change. A box that was large enough at the start may become too small after 100 ns if the protein extends.

The most common source of non-reproducible MD results. A system that hasn’t equilibrated is not representative of physiological conditions — temperature, pressure, and density are still drifting, and the protein may be in a metastable state far from its true equilibrium structure. Running production MD from a poorly equilibrated system means collecting data from an artifact rather than a real simulation.

The telltale sign: RMSD that rises steadily throughout what you’re calling the “production” phase. If RMSD hasn’t plateaued, the system hasn’t equilibrated — regardless of how long the equilibration phases nominally ran.

Particle Mesh Ewald (PME) — the standard long-range electrostatics method — requires the simulation box to be charge-neutral. A system with a non-zero net charge generates a divergent electrostatic energy under PBC that is handled by adding a uniform background charge, introducing a systematic error in the electrostatic calculation. The simulation runs and produces no error message, but all electrostatic energies are shifted by an unquantified constant.

This affects binding free energies, protein-ligand electrostatic interactions, and ion binding. It is especially problematic when the charged species is the ligand — a charged ligand in an uncorrected system will have spurious electrostatic interactions with the background charge.

Every biomolecular force field is parameterized against a specific water model. The protein parameters and water parameters are interdependent — they were fitted together to reproduce experimental observables including protein-water interactions, hydration free energies, and structural properties. Using a mismatched combination (CHARMM36 protein with TIP4P water, for example) introduces systematic errors in protein-water interactions that propagate into every property you measure.

Protein force fields cover the 20 amino acids. They say nothing about arbitrary drug molecules. Using protein-derived atom types for a ligand, mixing GAFF2 ligand parameters with CHARMM36 protein parameters, or using poorly validated ACPYPE output without checking the atom type assignments all produce physically wrong interactions between the ligand and protein. The simulation runs, but the binding affinity, ligand dynamics, and protein-ligand contacts are all artifacts.

The standard 2 femtosecond timestep is safe when X-H bonds are constrained with LINCS or SHAKE. Without these constraints, the fastest atomic motions (X-H bond vibrations, ~10 fs period) require a timestep of 1 fs or smaller. Using 2 fs without hydrogen constraints can cause numerical instability that either crashes the simulation with a Stepsize too small error or, more insidiously, produces a trajectory that looks normal but has accumulated integration errors.

The opposite problem also exists: some researchers use very long timesteps (4–5 fs with virtual sites) without properly validating that structural properties are preserved. A timestep that’s too long for the system’s fastest degrees of freedom produces systematic errors that don’t announce themselves.

The Berendsen thermostat and barostat equilibrate temperature and pressure quickly and stably — which is why they’re appropriate for equilibration. But they do not produce a correct thermodynamic ensemble: they suppress the natural fluctuations that define the canonical (NVT) or isothermal-isobaric (NPT) ensemble. Using Berendsen for production MD gives the wrong heat capacity, wrong compressibility, wrong free energies, and potentially wrong conformational populations. The simulation looks fine; the thermodynamics are wrong.

Periodic boundary conditions cause atoms to appear to teleport across the box boundary in the raw trajectory. This is not a physical artifact — it’s a coordinate wrapping convention. But it causes real problems when you run analysis on the raw trajectory: RMSD calculations give wrong values when atoms jump across boundaries, hydrogen bond distances appear impossibly large, and visualization shows a fragmented protein. The simulation is fine; the analysis is wrong.

In NVT simulations with no thermostat coupling (or very weak coupling), total energy should be approximately conserved over time. Systematic energy drift — the total energy rising or falling steadily throughout the simulation — indicates a problem with the integrator, timestep, or force field parameters. It’s a diagnostic signal that most researchers never check, leaving subtle simulation problems undetected.

Energy drift is measured in kJ/mol/atom/ns. A drift rate below 0.01 kJ/mol/atom/ns is excellent. Above 0.1 suggests a problem worth investigating.

A 100 ns simulation that samples a local energy minimum looks identical to a 100 ns simulation that has genuinely explored the relevant conformational space — unless you check convergence. Conclusions drawn from unconverged simulations are simply observations of where the system happened to be, not statements about equilibrium properties. This is the most common way MD studies overstate their confidence in results.

Common convergence checks include: comparing RMSF profiles from the first and second halves of the trajectory (they should be similar), checking whether RMSD has genuinely plateaued or is still slowly rising, and running multiple independent replicas with different initial velocities to confirm key observations are reproducible.