How Molecular Dynamics Simulation Works: Force Fields to Trajectories

A crystal structure or an AlphaFold model gives you one frozen pose of a protein, but proteins are not statues. They breathe, loops flap open and shut, side chains rotate, and a drug slides into a pocket that only exists for a few nanoseconds. Molecular dynamics (MD) is how you put that static coordinate file into motion and watch it wiggle under Newton's laws. The catch is that "put it into motion" hides a stack of decisions -- which force field, how to minimize, how to equilibrate temperature and pressure -- and getting any of them wrong produces a trajectory that looks fine in a movie but is physically meaningless. This post walks through what actually happens inside an MD run, the same way we walk through read alignment or variant calling.

The Core Idea: Newton's Laws on Atoms

At its heart MD is embarrassingly simple. Every atom is a point mass with a position and a velocity. You compute the force on each atom, use F = ma to get its acceleration, and step everything forward by a tiny slice of time. Repeat a few hundred million times and you have a trajectory: a movie of the system evolving.

The integration step is done with an algorithm like the leap-frog or velocity-Verlet integrator. The time step is the single most important number in the whole simulation. A typical value is 2 fs (femtoseconds). Why so small? Because the fastest motions in the system -- bonds involving hydrogen stretching -- vibrate on the order of 10 fs. To resolve a vibration you need several steps per period, so the time step has to stay well below the fastest motion or the integrator explodes. This is exactly why constraint algorithms like LINCS or SHAKE matter: they freeze the hydrogen bond lengths so you can safely use 2 fs instead of 0.5 fs and get four times the sampling for free.

Force Fields: Where the Physics Lives

The force in F = ma comes from the potential energy function, and that function is the force field. It is a set of parameters and equations describing how atoms interact: bond stretching, angle bending, dihedral torsions, plus the non-bonded terms -- van der Waals (Lennard-Jones) and electrostatics (Coulomb). A force field is not derived from first principles; it is fit to quantum calculations and experimental data. Choosing one is a scientific decision, not a default to accept blindly.

For proteins in 2025-2026, two families dominate:

The pairing column matters more than beginners expect. ff19SB was parameterized against the OPC water model, and using it with plain TIP3P degrades the very backbone improvements you chose it for. If your system has a small-molecule ligand, the protein force field does not cover it -- you generate ligand parameters separately with something like GAFF2 and antechamber, or the CHARMM General Force Field (CGenFF).

Building the System: Solvation and Ions

Before any dynamics run, you build a simulation box. A protein floating in vacuum is unphysical, so you place it in a box of explicit water molecules -- this is where the water model choice becomes concrete. You leave a margin (often 1.0 nm to 1.2 nm) between the protein and the box edge so the molecule never "sees" its own periodic image across the boundary.

Real cells are not distilled water, so you add ions to neutralize the system's net charge and often to reach a physiological concentration around 0.15 M NaCl. In GROMACS this is the gmx solvate and gmx genion step; in AMBER's tleap it is solvateBox followed by addIons. Periodic boundary conditions then make the box tile infinitely, so an atom leaving the right face re-enters on the left, eliminating hard walls.

Energy Minimization: Removing the Clashes

Freshly built systems are full of tiny geometric horrors -- two atoms placed a fraction of an angstrom too close, a water oxygen overlapping a side chain. If you launched dynamics immediately, those clashes would produce enormous forces, the integrator would blow up, and you would see the dreaded LINCS warning or a segmentation fault.

Energy minimization fixes this before adding any kinetic energy. It is not dynamics at all; it is optimization. Algorithms like steepest descent walk the system downhill on the potential energy surface, nudging atoms until the maximum force drops below a threshold such as 1000 kJ/mol/nm. Minimization does not find the global minimum and is not supposed to -- it just removes the worst clashes so dynamics can start safely. Think of it as relaxing the geometry, not folding the protein.

Equilibration: NVT then NPT

Now you add motion, but carefully, in two stages. This is the step beginners skip and later regret.

In the NVT ensemble you hold the number of particles, volume, and temperature fixed. You assign initial velocities from a Maxwell-Boltzmann distribution at your target temperature, typically 300 K, and run a short simulation (often 100 ps to 500 ps) with a thermostat -- the V-rescale (Bussi) or Nosé-Hoover thermostat -- that gently couples the system to a heat bath until the temperature stabilizes. During NVT you usually keep position restraints on the protein heavy atoms so the solvent equilibrates around a fixed solute.

In the NPT ensemble you swap volume for pressure, adding a barostat -- Parrinello-Rahman or the C-rescale barostat -- to hold pressure near 1 bar. This lets the box volume adjust so the solvent reaches the correct density (around 1 g/cm³ for water). You gradually release the position restraints. Only when both temperature and pressure have plateaued, and density has converged, is the system truly equilibrated.

Production Dynamics: The Trajectory You Analyze

With restraints gone and the ensemble stable, you launch the production run -- the long, unrestrained simulation that actually samples the physics you care about. Lengths depend on the question: 100 ns for a quick look at side-chain flexibility, hundreds of nanoseconds to microseconds for conformational changes or ligand residence. The engine writes coordinates to a trajectory file every few picoseconds (GROMACS .xtc, AMBER .nc), and that trajectory is your raw data, the MD equivalent of a BAM file.

Analysis then extracts signal: RMSD to check structural stability, RMSF per residue to map flexibility, radius of gyration for compactness, hydrogen-bond occupancy, and free-energy estimates via MM-PBSA or MM-GBSA for binding. A flat RMSD plateau suggests a stable fold; a still-climbing RMSD at the end of your run means you have not sampled long enough.

Choosing an Engine: GROMACS, AMBER, or OpenMM

The three engines you will meet most often solve the same equations but feel very different in practice.

GROMACS is the usual starting point because it is free, fast, and has the deepest tutorial ecosystem. OpenMM wins when you want to script the whole workflow in Python or plug in a machine-learning potential. AMBER is favored when you specifically want its native force fields validated end to end. Note that force fields and engines are largely portable -- you can run CHARMM36m in GROMACS or ff19SB in OpenMM -- so the engine choice is more about workflow and hardware than physics.

If you have not set up a compute environment yet, MD tools install cleanly through conda; our guides on conda environments and Bioconda channels and installing Miniconda on Apple Silicon cover the foundation.

Practical Takeaways

• MD is Newton's laws on atoms: pick a small time step (2 fs with LINCS or SHAKE constraints on hydrogens) or the integrator will explode.
• The force field is a scientific choice, not a default. ff19SB pairs with OPC water; CHARMM36m pairs with modified TIP3P. Never mix families.
• Always run the full pipeline in order: build and solvate, minimize to kill clashes, equilibrate under NVT then NPT, then production. Skipping equilibration corrupts everything downstream.
• Check that temperature and density have plateaued before trusting a production trajectory, and check that RMSD has stabilized before drawing conclusions.
• Treat the trajectory like sequencing data: it is raw output that needs analysis (RMSD, RMSF, MM-PBSA), not an answer by itself.
• Start with GROMACS for tutorials, reach for OpenMM when you want Python control, and use AMBER for its native force fields.

Molecular dynamics rewards patience and punishes shortcuts, but the payoff is a physically grounded view of how a protein actually moves that no static structure can give you. Once the concepts here feel solid, the natural next step is a hands-on run with a real ligand -- see our GROMACS protein-ligand simulation tutorial, and pair it with molecular docking in AutoDock Vina to generate the starting poses your simulation refines.