Written and maintained by CASRAI Editorial Board
Published
On this page
Short answer
Molecular dynamics simulation is a computational method that calculates the time-dependent physical trajectories of atoms and molecules by numerically integrating classical equations of motion under an empirical force field. By computing interatomic forces at femtosecond increments, simulations reveal conformational dynamics, thermodynamic properties, and biomolecular interaction mechanisms that remain inaccessible to static structural biology or macroscopic laboratory experiments.
Definition and scope
Molecular dynamics (MD) simulation is a foundational computational technique within physical chemistry, structural biophysics, and materials science. Unlike static structural determination methods—such as X-ray crystallography, cryo-electron microscopy, or predictive modeling from the AlphaFold protein structure database—molecular dynamics models the continuous conformational landscape of molecular systems across temporal regimes ranging from femtoseconds to microseconds.
The operational scope spans all-atom models, where every heavy atom, hydrogen, and explicit solvent molecule is represented, to coarse-grained representations that group atoms into interaction beads to capture large macromolecular assemblies and extended timescales.
How it works or how it is used in practice
Classical molecular dynamics replaces quantum mechanical wavefunctions with empirical potential energy functions known as force fields. The total potential energy U(r) is partitioned into bonded and non-bonded terms:
U(r) = U_bonded + U_nonbonded
Bonded interactions maintain covalent structure via harmonic bond stretching (U_bond = 0.5 × k_b × (r – r_0)²), harmonic angle bending (U_angle = 0.5 × k_θ × (θ – θ_0)²), and periodic dihedral torsions (U_torsion = 0.5 × V_n × [1 + cos(nϕ – γ)]). Non-bonded interactions include Lennard-Jones 12-6 potentials (U_LJ = 4ε × [(σ / r_ij)¹² – (σ / r_ij)⁶]) for Pauli repulsion and van der Waals dispersion, alongside Coulomb electrostatics (U_coulomb = (q_i × q_j) / (4πε_0 × r_ij)) evaluated using Particle Mesh Ewald (PME) summation in periodic systems.
Numerical Integration and Time Steps
Atomic trajectories are generated by integrating Newton’s second law, F_i = m_i × a_i = -∇_i U(r), using algorithms such as Leap-Frog or Velocity Verlet. Because high-frequency covalent bonds involving hydrogen vibrate on a 10-femtosecond timescale, numerical stability dictates an unconstrained integration time step of 1 femtosecond (fs) or less. Applying rigid constraint algorithms (such as SHAKE or LINCS) fixes bond lengths to hydrogen, safely extending the integration time step to 2 fs. Hydrogen mass repartitioning (HMR) rebalances atomic masses to permit production time steps up to 4 fs.
Standard Simulation Execution Protocol
A standard simulation workflow follows six sequential steps:
- 1. System Preparation: Assemble atomic coordinates, resolve missing loops, and assign residue protonation states matching target pH.
- 2. Solvation and Neutralization: Enclose the molecule in a periodic solvent box with a minimum 1.0 nm buffer, solvate with explicit water (e.g., TIP3P), and add counter-ions (Na⁺, Cl⁻) to achieve electrical neutrality and physiological ionic strength (~0.15 M NaCl).
- 3. Energy Minimization: Run Steepest Descents minimization until maximum forces drop below 1,000 kJ/(mol × nm) to eliminate steric clashes.
- 4. NVT Thermalization: Heat the system to target temperature (e.g., 300 K or 310 K) under constant volume using velocity rescaling or Langevin thermostats with position restraints on solute heavy atoms.
- 5. NPT Pressurization: Equilibrate system density and box volume at 1.0 bar using a barostat (such as Parrinello-Rahman) while gradually releasing solute restraints.
- 6. Production Run: Propagate the unconstrained trajectory over planned durations, sampling coordinates at defined intervals.
Types or related concepts
Simulations are defined by their statistical thermodynamic ensemble. The table below outlines standard ensembles, control mechanisms, and research applications:
| Ensemble Name | Fixed State Variables | Coupling Mechanism | Primary Research Application | Standard Reference |
|---|---|---|---|---|
| Microcanonical (NVE) | Particle number (N), Volume (V), Total energy (E) | Isolated system (no thermostat or barostat) | Validating numerical integrator stability and verifying energy conservation without drift. | Braun et al., 2019 |
| Canonical (NVT) | Particle number (N), Volume (V), Temperature (T) | Nosé-Hoover, Langevin, or Bussi-Donadio-Parrinello thermostat | Initial thermal equilibration, constant-volume solvent systems, and fixed-area membrane interfaces. | Hollingsworth and Dror, 2018 |
| Isothermal-Isobaric (NPT) | Particle number (N), Pressure (P), Temperature (T) | Thermostat coupled with Parrinello-Rahman or Monte Carlo barostat | Standard biomolecular production simulations, solvated ligand binding, and liquid-state density evaluations. | Braun et al., 2019 |
| Grand Canonical (μVT) | Chemical potential (μ), Volume (V), Temperature (T) | Particle insertion/deletion schemes coupled to thermal bath | Porous material adsorption, ion channel permeation, and chemical equilibrium studies. | GROMACS Manual, 2024 |
Molecular dynamics interfaces with complementary computational methods:
- Monte Carlo Simulation: Unlike molecular dynamics, which computes deterministic, time-ordered pathways from forces, Monte Carlo simulation samples equilibrium states through stochastic probability evaluations without physical time tracking.
- Molecular Docking: While molecular docking software predicts static ligand binding poses, molecular dynamics captures receptor flexibility, solvent effects, and induced-fit binding stability.
Where it matters for labs, clinics or supply decisions
Computational simulations guide experimental design, screening, and high-performance resource allocation across academic and industrial laboratories.
GPU Hardware Allocation
Modern simulation packages (GROMACS, AMBER, OpenMM) leverage GPU acceleration. Because throughput is governed primarily by single-precision floating-point throughput and memory bandwidth, modern workstation GPUs routinely deliver 150 to 300 nanoseconds per day for typical 50,000-atom solvated systems, enabling rapid in silico screening prior to experimental synthesis.
Data Stewardship and Trajectory Storage
Production trajectories generate massive datasets. Storing coordinate trajectories in compressed binary formats (such as XTC) with 0.001 nm precision reduces storage footprints by 70% to 90% compared to uncompressed raw files. In accordance with the RDA FAIR Data Maturity Model, reproducible data publishing requires depositing starting coordinates, force field parameters, and run scripts in open repositories.
Common misunderstandings
- Misunderstanding: Extended simulation length compensates for poor starting structures. Unresolved steric clashes, incorrect protonation states, or erroneous disulfide bonds persist and propagate throughout dynamics. Initial structural validation remains indispensable.
- Misunderstanding: Classical force fields simulate covalent bond breaking. Standard force fields employ harmonic bonded potentials that cannot model chemical bond rupture or formation. Modeling covalent reactions requires QM/MM or reactive force fields.
- Misunderstanding: Truncating electrostatic interactions with simple spherical cutoffs is harmless. Abrupt cutoffs generate artificial forces at boundary interfaces, introducing structural artifacts and spurious heating. Long-range grid electrostatics (PME) are necessary.
- Misunderstanding: The Berendsen thermostat samples a true canonical ensemble. While effective for rapid initial thermal relaxation, Berendsen velocity scaling suppresses physical kinetic energy fluctuations. True NVT sampling requires Langevin dynamics or Nosé-Hoover thermostats.
- Misunderstanding: Unconstrained systems can safely use 2-femtosecond time steps. High-frequency covalent vibrations involving hydrogen destabilize unconstrained integration above 1 fs. Bond-length constraint algorithms (LINCS or SHAKE) are mandatory for 2 fs time steps.
FAQ
What is the difference between molecular dynamics and Monte Carlo simulation?
Molecular dynamics calculates deterministic, time-resolved trajectories by integrating Newton’s equations of motion from interatomic forces, capturing transport kinetics and physical pathways. Monte Carlo simulation samples equilibrium configurations through stochastic state moves based on the Boltzmann factor, evaluating thermodynamic averages without time progression.
Why are algorithms like LINCS or SHAKE necessary for 2-femtosecond time steps?
Covalent bonds to hydrogen vibrate at approximately 10 femtoseconds. To maintain numerical integration stability, time steps must remain below roughly one-tenth of the fastest vibrational period (≤ 1 fs). Algorithms like LINCS and SHAKE constrain these high-frequency bond lengths, eliminating high-frequency oscillations and permitting stable 2 fs time steps.
How do researchers decide between NVT and NPT ensembles?
The NVT ensemble is selected for initial heating and systems where volume or interfacial surface area must stay fixed, such as surface adsorption or monolayer studies. The NPT ensemble is used for production runs of solvated macromolecules because allowing the box volume to adjust maintains natural atmospheric pressure and physiological solvent density.
What is the difference between all-atom and coarse-grained molecular dynamics?
All-atom simulations model every individual atom, providing atomic-level chemical detail for hydrogen bonding and ligand contacts over sub-microsecond timescales. Coarse-grained simulations group clusters of atoms into single interaction beads, reducing computational degrees of freedom to access multi-microsecond or millisecond biological processes.
Sources
- Braun, E., Moosavi, S. M., and Smit, B. (2019). Best Practices for Foundations in Molecular Simulations [Article v1.0]. Living Journal of Computational Molecular Science, 1(1), 5957. doi:10.33011/livecoms.1.1.5957. PMCID: PMC6884151. Read November 11, 2026.
- Hollingsworth, S. A., and Dror, R. O. (2018). Molecular Dynamics Simulation for All. Neuron, 99(6), 1129–1143. doi:10.1016/j.neuron.2018.08.011. PMCID: PMC6209097. Read November 11, 2026.
- GROMACS Development Team. (2024). GROMACS Reference Manual. Algorithms and Molecular Dynamics Integrators. URL: https://manual.gromacs.org/current/reference-manual/algorithms/molecular-dynamics.html. Read November 11, 2026.
