Molecular Dynamics: Equations of Motion

Introduction

Molecular dynamics, or MD, turns a potential-energy function into a time-dependent trajectory. Given positions and velocities, the force field supplies forces; a numerical integrator advances the system by a short time step. Repeating this calculation can show atomic fluctuations, diffusion and conformational transitions. It is tempting to treat the trajectory as a film of exact molecular motion, but its path depends on the force field, chosen time step, initial conditions and boundary treatment. Understanding the equations and integration error helps separate numerical reliability from chemical validity.

Core explanation

For atom i with mass mi, Newton's equation is mi d²ri/dt² = Fi. The force is Fi = −∇i U(r₁, r₂, …), where U is the selected potential energy. At each step the program evaluates the current forces, updates velocities and positions, and repeats. In a constant-energy idealization with a time-independent conservative potential and no thermostat, the total energy E = kinetic + potential should remain approximately conserved. Finite steps introduce integration error; force cutoffs, constraints and numerical precision can also matter. The simulation is classical for nuclear motion and commonly uses a classical force field, though the forces can in principle come from electronic-structure calculations in ab initio MD.

A basic Verlet-family integrator is favored because it is stable and approximately time reversible for suitable steps. In leap-frog form, velocities live at half steps: v(t + Δt/2) = v(t − Δt/2) + F(t)Δt/m, then r(t + Δt) = r(t) + v(t + Δt/2)Δt. Velocity Verlet instead stores positions and velocities together at whole steps while using half-step velocity updates. The OpenMM integrator documentation specifies the leap-frog equations and warns that position and velocity time alignment matters when evaluating energy. These algorithms approximate continuous motion; they do not make an inadequate force field more physically correct.

The time step must resolve the fastest relevant motions. X–H bond stretches are rapid, so many biomolecular simulations constrain selected bonds involving hydrogen and use a somewhat larger step than an unconstrained model. There is no universal safe step: stiff forces, high temperatures, reactive collisions or a different integrator can require a smaller one. A trajectory that does not visibly explode can still have biased energy, pressure or dynamics. Test a shorter step and inspect energy conservation in a suitable constant-energy diagnostic. GROMACS's molecular-dynamics algorithm manual describes integration and the role of constraints and coupling.

Initial velocities are often sampled from a temperature-related distribution, but one set of initial velocities does not uniquely represent the ensemble. An equilibration period allows the system to relax from constructed starting conditions. Production sampling should follow stable temperature, density and energy behavior appropriate to the chosen ensemble. Time averages approximate equilibrium expectations only if sampling is adequate and the system explores the relevant states. A long trajectory trapped in one conformation can be less informative about equilibrium than shorter independent trajectories that visit all relevant basins.

Forces need to be consistent with the potential-energy function. If a software implementation reports U but applies a mismatched gradient, an NVE energy drift test can reveal trouble. Long-range electrostatics, cutoff schemes and periodic boundaries affect forces and must be chosen consistently. A thermostat or barostat intentionally exchanges energy with an external mathematical reservoir, so total physical-system energy need not be constant during such runs; use the diagnostic appropriate to that ensemble. A stable NVE energy trace supports numerical integration, but it does not validate the chemistry of the potential.

Step-by-step reasoning

1. Define the atoms, masses, potential-energy model, topology and boundary conditions. 2. Minimize severe initial overlaps before assigning sensible initial velocities. 3. Choose an integrator, time step and any bond constraints based on the fastest motions. 4. Run a short diagnostic with a smaller step and, where appropriate, check constant-energy behavior. 5. Equilibrate under the intended thermodynamic conditions, then collect production trajectories. 6. Analyze independent samples and distinguish numerical convergence from force-field and sampling error.

Visual explanation

Draw a time axis marked t, t + Δt/2 and t + Δt. Put position dots at whole steps and velocity arrows at half steps to illustrate leap-frog integration. Above it, draw a ball moving on a curved potential-energy surface; the local slope gives force. A coarse-step path visibly cuts across the curvature, while a finer-step path follows it more closely. Beneath, show a nearly flat total-energy trace for a well-behaved NVE test and a drifting trace that calls for investigation.

Real-world analogy

Following a winding hiking path by taking short compass steps traces it more accurately than taking long straight jumps between occasional readings. An MD integrator samples forces at discrete times and makes short updates. However, even tiny steps cannot fix an incorrect map of the terrain: integration accuracy and force-field accuracy are separate questions.

Real-world example

A researcher simulates a small protein in water. An initial structure contains overlapping solvent molecules, so the system is minimized and equilibrated before data collection. With constraints on suitable hydrogen bonds, the researcher tests the selected time step against a shorter one, checks stable ensemble properties and launches independent seeds. One run remains folded while another visits an alternate loop conformation. Reporting only the first trajectory as the protein's unique behavior would confuse initial-condition dependence with equilibrium evidence.

Why?

Why can a too-large time step create energy drift? A numerical update assumes forces and trajectories can be approximated over Δt. If the step skips over rapid curvature changes, the discrete path no longer closely follows the continuous equations of motion. Errors can accumulate or destabilize the trajectory. Shortening Δt generally improves integration but raises cost because more force evaluations are needed for the same physical duration.

Common misconception

“A longer trajectory is always more accurate.” It may remain trapped or use a poor force field. “A thermostat fixes an unstable time step.” Temperature control can conceal but not cure bad integration. “MD trajectories are exact quantum paths of atoms.” Ordinary MD uses classical nuclear dynamics and an approximate potential. “Constant energy in one test proves the force field is chemically correct.” It checks numerical consistency, not agreement with experiment.

Worked example

Suppose a simulation advances 5,000,000 steps at 2 femtoseconds per step. The total simulated time is 10,000,000 fs = 10 ns. Reducing the step to 1 fs while keeping 5,000,000 steps gives only 5 ns, although integration may be more accurate. To simulate the same 10 ns at 1 fs requires 10,000,000 steps and roughly twice as many force evaluations. If an NVE diagnostic shows a steady 5-kJ mol⁻¹ energy rise over the first nanosecond at 2 fs but a nearly flat trace at 1 fs, the larger step is suspect for that model. The numerical thresholds depend on system size and method; these illustrative values are not universal tolerances.

Quick check

1. What is the force on atom i in an energy-based molecular-mechanics model? Answer: The negative gradient of potential energy with respect to that atom's coordinates. 2. Does an NVT simulation have to conserve the simulated system's kinetic-plus-potential energy? Answer: No. A thermostat exchanges energy to maintain the intended temperature distribution.

Exam focus

Write mi ai = Fi and Fi = −∇i U, then explain how discrete integration advances a trajectory. Convert steps and femtoseconds into nanoseconds correctly. State why the fastest motions limit Δt and why constraints can change the practical choice. Distinguish energy-conservation diagnostics in NVE from temperature control in NVT. Explain why numerical stability, adequate sampling and force-field validity are separate requirements.

Advanced insight

Symplectic, reversible integrators can show bounded oscillations around a nearby conserved quantity rather than perfect conservation of the continuous Hamiltonian. Floating-point precision, constraint tolerances and long-range-force update schedules can still influence measured drift. Stochastic integrators deliberately add random and friction forces, so their pathwise energy behavior differs from deterministic NVE Verlet. For kinetic interpretation, thermostat parameters and friction can alter time correlations even when equilibrium distributions are approximately correct. The integration scheme should be chosen for the observable, not only for maximum step size.

Summary

MD propagates positions and velocities under forces derived from a chosen potential. Verlet-type finite-step algorithms approximate Newtonian motion efficiently, but time step and constraint choices must resolve relevant dynamics. Energy and stability tests detect numerical problems; equilibration and independent sampling address ensemble quality. A well-integrated trajectory remains limited by its force field and physical model, so conclusions require all three checks: numerical behavior, sampling and chemical validation.

Practice questions

1. How long is a run of 2,000,000 steps at 1 fs each? Answer: 2,000,000 fs = 2 ns. 2. If the potential increases in the positive x direction, which way does the x force point? Answer: Toward negative x, since F x = −∂U/∂x. 3. Why may constraining X–H bond vibrations permit a longer step? Answer: The fastest bond-stretch motion is removed from explicit integration, reducing the highest frequency that must be resolved. 4. Does a flat NVE energy trace establish that a ligand binds correctly to a protein? Answer: No. It supports integration consistency; binding accuracy depends on force-field quality and sampling.