Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Introduction to Molecular Dynamics: Concepts, Algorithms, and Formalism

1. The Philosophy: A Computational Microscope

1.1 Two Paradigms of Molecular Simulation

In previous sections, we used Statistical Mechanics to calculate properties based on the probability of states. We asked: “In an infinite ensemble of copies, what is the average energy?”

Molecular Dynamics (MD) takes a fundamentally different approach. It asks: “How does this specific system evolve over time?”

MD is a deterministic simulation of the classical Many-Body problem. By solving Newton’s laws of motion for every atom, we generate a Trajectory: a chronological movie of the microscopic world.

ApproachQuestion AskedOutputTime Information
Statistical MechanicsWhat are equilibrium properties?Ensemble averages (T, P, G, S)No
Molecular DynamicsHow does the system evolve?Trajectories, mechanismsYes

1.2 What Can MD Tell Us?

MD simulations provide unique insights that cannot be obtained from experiments or static calculations:

  1. Dynamical Properties:

    • Diffusion coefficients

    • Viscosity and thermal conductivity

    • Vibrational spectra

    • Reaction rates and transition states

  2. Microscopic Mechanisms:

    • Protein folding pathways

    • Crystal nucleation

    • Crack propagation in materials

    • Drug binding mechanisms

  3. Structural Fluctuations:

    • Conformational changes

    • Phase transitions

    • Defect formation and migration

1.3 Limitations and Approximations

MD simulations, while powerful, have important limitations:

2. The Formalism: Equations of Motion

2.1 Phase Space Representation

Formally, the state of a system with NN atoms is defined by a single point in Phase Space (Γ\Gamma) with 6N6N dimensions:

The trajectory of the system traces out a path in this high-dimensional space. Each point uniquely defines the microstate of the system.

Key Insight: In classical mechanics, if you know the exact positions and momenta at time tt, you can predict them at all future times (determinism).

2.2 The Hamiltonian Formulation

The evolution of this point is governed by the Hamiltonian (H\mathcal{H}), which represents the total energy of the system.

H(rN,pN)=K(pN)+U(rN)\mathcal{H}(\mathbf{r}^N, \mathbf{p}^N) = \mathcal{K}(\mathbf{p}^N) + \mathcal{U}(\mathbf{r}^N)
(1)

Where:

2.3 Hamilton’s Equations of Motion

The motion is derived from Hamilton’s equations:

  1. Velocity: r˙i=∂H∂pi=pimi\dot{\mathbf{r}}_i = \frac{\partial \mathcal{H}}{\partial \mathbf{p}_i} = \frac{\mathbf{p}_i}{m_i}

  2. Force: p˙i=−∂H∂ri=Fi\dot{\mathbf{p}}_i = -\frac{\partial \mathcal{H}}{\partial \mathbf{r}_i} = \mathbf{F}_i

These are equivalent to Newton’s second law: Fi=miai=mir¨i\mathbf{F}_i = m_i \mathbf{a}_i = m_i \ddot{\mathbf{r}}_i

Conceptual Check: The Conservative Force

The second equation is critical. It tells us that Force is the negative gradient of Potential Energy.

Fi=−∇iU(r1,...,rN)=−∂U∂ri\mathbf{F}_i = -\nabla_i U(\mathbf{r}_1, ..., \mathbf{r}_N) = -\frac{\partial U}{\partial \mathbf{r}_i}
(2)

This means atoms always “roll downhill” on the Potential Energy Surface, trying to minimize their potential energy. The steeper the slope, the larger the force.

2.4 Conservation Laws

Hamilton’s equations automatically conserve several important quantities:

These conservation laws are critical for the stability and physical validity of MD simulations.

3. The Algorithm: Anatomy of an MD Code

3.1 The Basic Structure

How do we translate this formalism into code? Despite differences between software like LAMMPS, GROMACS, NAMD, or AMBER, they all share a common algorithmic “skeleton.”

The simulation proceeds in a loop that advances time in small, discrete steps (Δt\Delta t), typically 1-2 femtoseconds (fs).

3.2 Pseudo-code Representation

# ========================================
# Initialize the System
# ========================================
t = 0
Initialize Positions (x) from structure file or lattice
Initialize Velocities (v) from Maxwell-Boltzmann distribution
Initialize Forces (F) by calling force calculation

# ========================================
# The Main MD Loop
# ========================================
while (t < t_max):
    # 1. Calculate Forces (The "Expensive" Step)
    #    This takes ~90% of computational time
    forces = calculate_forces(x)
    
    # 2. Integrate Equations of Motion (Verlet)
    #    Update positions and velocities
    x_new, v_new = velocity_verlet(x, v, forces, dt)
    
    # 3. Apply Constraints (if any)
    #    E.g., SHAKE/RATTLE for bond lengths
    apply_constraints(x_new, v_new)
    
    # 4. Apply Thermostat/Barostat (if needed)
    #    Control temperature and/or pressure
    if (using_thermostat):
        v_new = thermostat_update(v_new, T_target)
    if (using_barostat):
        x_new, box_size = barostat_update(x_new, P_target)
    
    # 5. Sampling (Thermodynamics)
    #    Calculate and store properties
    if (t % sample_freq == 0):
        calculate_properties(x_new, v_new)
        write_trajectory(x_new, v_new)
        
    # 6. Advance Time
    t = t + dt
    x, v = x_new, v_new
end loop

# ========================================
# Finalize
# ========================================
calculate_averages()
write_final_configuration()

3.3 Component 1: Initialization

We need to set the starting point in Phase Space.

Positions:

Velocities:

Initial Force Calculation:

3.4 Timestep Selection

The choice of timestep Δt\Delta t is critical:

Rule of Thumb: Δt\Delta t should be ~1/10 of the fastest vibrational period in the system

System TypeTypical Timestep
All-atom with hydrogens1 fs
All-atom, constrained H bonds2 fs
United atom (CH₃ groups)2-5 fs
Coarse-grained10-50 fs

4. The Force Problem: Interatomic Potentials

4.1 The Central Challenge

The accuracy of any MD simulation depends entirely on the quality of the Potential Energy Surface, U(rN)\mathcal{U}(\mathbf{r}^N). Since solving the Schrödinger equation for every step is too expensive (quantum chemistry scales as N3N^3 to N7N^7), we use Empirical Force Fields.

4.2 The Car Engine Analogy

Imagine a car driving on a straight road. Newton’s laws tell us how it moves (kinematics), but to calculate its acceleration, we need to know the energy produced by the engine (dynamics).

Similarly, the Interatomic Potential is the driving engine of the simulation. The integrator (Verlet) is the transmission, but without accurate forces, the trajectory is meaningless.

4.3 The Lennard-Jones (12-6) Potential

For noble gases and simple fluids, the standard model is the Lennard-Jones potential:

ULJ(r)=4ϵ[(σr)12−(σr)6]U_{LJ}(r) = 4\epsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right]
(4)

Where:

Physical Interpretation:

The Force:

FLJ(r)=−dUdr=24ϵσ[2(σr)13−(σr)7]F_{LJ}(r) = -\frac{dU}{dr} = \frac{24\epsilon}{\sigma}\left[2\left(\frac{\sigma}{r}\right)^{13} - \left(\frac{\sigma}{r}\right)^{7}\right]
(5)

4.4 Python: Visualizing the Physics

import numpy as np
import matplotlib.pyplot as plt

# LJ Parameters for Argon
epsilon = 1.0  # Normalized energy (119.8 K in real units)
sigma = 1.0    # Normalized distance (3.405 Å in real units)

r = np.linspace(0.9, 3.0, 200)

# Potential U(r)
U = 4 * epsilon * ((sigma/r)**12 - (sigma/r)**6)

# Force F(r) = -dU/dr
F = 24 * epsilon / sigma * (2*(sigma/r)**13 - (sigma/r)**7)

# Create figure with two subplots
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 10))

# Plot potential
ax1.plot(r, U, 'b-', linewidth=2, label='U(r)')
ax1.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax1.axhline(y=-epsilon, color='r', linestyle='--', alpha=0.5, label=r'$-\epsilon$')
ax1.axvline(x=2**(1/6)*sigma, color='g', linestyle='--', alpha=0.5, label=r'$r_{min}$')
ax1.set_ylim(-1.5, 2)
ax1.set_xlabel('Distance r/σ', fontsize=12)
ax1.set_ylabel('Energy U/ε', fontsize=12)
ax1.set_title("Lennard-Jones Potential", fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(alpha=0.3)

# Plot force
ax2.plot(r, F, 'r-', linewidth=2, label='F(r)')
ax2.axhline(y=0, color='k', linestyle='--', alpha=0.3)
ax2.axvline(x=2**(1/6)*sigma, color='g', linestyle='--', alpha=0.5, label=r'$r_{min}$')
ax2.set_xlabel('Distance r/σ', fontsize=12)
ax2.set_ylabel('Force F×σ/ε', fontsize=12)
ax2.set_title("Lennard-Jones Force", fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('LJ_potential_force.png', dpi=150)
plt.show()

4.5 Cutoffs and Long-Range Corrections

Computing all pairwise interactions scales as O(N2)O(N^2). For N=106N = 10^6 atoms, this means 1012 evaluations per timestep!

Cutoff Radius (rcr_c):

Tail Correction: For LJ, we can analytically correct for neglected long-range interactions:

Utail=8πNρϵσ33[13(σrc)9−(σrc)3]U_{tail} = \frac{8\pi N \rho \epsilon \sigma^3}{3}\left[\frac{1}{3}\left(\frac{\sigma}{r_c}\right)^9 - \left(\frac{\sigma}{r_c}\right)^3\right]
(6)

Electrostatics: Coulomb interactions are long-range (∝1/r\propto 1/r) and cannot be simply cut off. We use:

5. Integration: Moving the Atoms

5.1 Why Not Simple Integration?

We have the initial state {r(t),v(t)}\{\mathbf{r}(t), \mathbf{v}(t)\} and the forces F(t)\mathbf{F}(t). Why can’t we use simple calculus?

Euler Method:

r(t+Δt)=r(t)+v(t)Δt\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t) \Delta t
(7)

v(t+Δt)=v(t)+F(t)mΔt\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{\mathbf{F}(t)}{m} \Delta t
(8)

Problems:

5.2 Required Properties of MD Integrators

A good integrator must be:

  1. Time-reversible: r(t+Δt)→r(t)\mathbf{r}(t + \Delta t) \rightarrow \mathbf{r}(t) when velocities are reversed

  2. Symplectic: Preserves phase space volume (Liouville’s theorem)

  3. Energy-conserving: Total energy fluctuates but doesn’t drift

  4. Efficient: Requires few force evaluations (forces are expensive!)

  5. Stable: Works for reasonably large timesteps

5.3 The Velocity Verlet Algorithm

The industry standard is Velocity Verlet, which satisfies all the above requirements.

Algorithm:

  1. Half-Kick: Update velocities by half a timestep using current forces.

    v(t+Δt2)=v(t)+F(t)2mΔt\mathbf{v}\left(t + \frac{\Delta t}{2}\right) = \mathbf{v}(t) + \frac{\mathbf{F}(t)}{2m}\Delta t
    (9)

  2. Drift: Update positions using the half-step velocities.

    r(t+Δt)=r(t)+v(t+Δt2)Δt\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}\left(t + \frac{\Delta t}{2}\right)\Delta t
    (10)

  3. Force Update: Calculate new forces F(t+Δt)\mathbf{F}(t + \Delta t) at the new positions r(t+Δt)\mathbf{r}(t + \Delta t).

  4. Half-Kick: Complete the velocity update.

    v(t+Δt)=v(t+Δt2)+F(t+Δt)2mΔt\mathbf{v}(t + \Delta t) = \mathbf{v}\left(t + \frac{\Delta t}{2}\right) + \frac{\mathbf{F}(t + \Delta t)}{2m}\Delta t
    (11)

Properties:

5.4 Python Implementation

import numpy as np

def velocity_verlet(positions, velocities, forces, masses, dt, force_function):
    """
    Velocity Verlet integrator
    
    Parameters:
    -----------
    positions : ndarray, shape (N, 3)
        Current positions
    velocities : ndarray, shape (N, 3)
        Current velocities
    forces : ndarray, shape (N, 3)
        Current forces
    masses : ndarray, shape (N,)
        Particle masses
    dt : float
        Timestep
    force_function : callable
        Function that computes forces from positions
        
    Returns:
    --------
    positions_new : ndarray
        Updated positions
    velocities_new : ndarray
        Updated velocities
    forces_new : ndarray
        Updated forces
    """
    # Half-step velocity update
    velocities_half = velocities + 0.5 * forces / masses[:, np.newaxis] * dt
    
    # Full-step position update
    positions_new = positions + velocities_half * dt
    
    # Calculate new forces
    forces_new = force_function(positions_new)
    
    # Half-step velocity update to complete the step
    velocities_new = velocities_half + 0.5 * forces_new / masses[:, np.newaxis] * dt
    
    return positions_new, velocities_new, forces_new

5.5 Alternative Integrators

Leapfrog Algorithm:

Runge-Kutta Methods:

Multiple Timestep (RESPA):

6. Controlling the Environment: Thermostats & Barostats

6.1 The Ensemble Problem

Newton’s laws naturally conserve Total Energy (E=K+U=constantE = K + U = \text{constant}). A standard MD simulation therefore generates the Microcanonical (NVE) Ensemble:

However, real experiments are usually done at:

Why do we care?

6.2 The Concept of Coupling

Imagine your simulation box is a test tube sitting in a large water bath (the reservoir). Heat flows in and out to keep TT constant. Similarly, a piston can compress/expand the box to keep PP constant.

Extended System Method: The reservoir itself becomes a dynamical variable with its own equation of motion. The combined system (simulation + reservoir) is isolated and conserves energy, but the simulation part exchanges energy with the reservoir.

6.3 Thermostats (NVT)

Temperature is related to kinetic energy via the equipartition theorem:

⟨K⟩=32NkBT\langle K \rangle = \frac{3}{2}Nk_B T
(12)

To control TT, we must manipulate the velocities.

6.3.1 Velocity Rescaling (Berendsen Thermostat)

Method: At each step, multiply all velocities by a scaling factor:

λ=1+ΔtτT(T0T(t)−1)\lambda = \sqrt{1 + \frac{\Delta t}{\tau_T}\left(\frac{T_0}{T(t)} - 1\right)}
(13)

Where:

Implementation:

vi(t+Δt)=λvi(t+Δt)\mathbf{v}_i(t + \Delta t) = \lambda \mathbf{v}_i(t + \Delta t)
(14)

Pros:

Cons:

6.3.2 Nosé-Hoover Thermostat

Method: Introduce an additional degree of freedom ss (the “heat bath”) with its own mass QQ and coordinate η\eta. The extended Hamiltonian is:

Hext=∑ipi22mi+U(rN)+pη22Q+NfkBT0η\mathcal{H}_{ext} = \sum_i \frac{p_i^2}{2m_i} + U(\mathbf{r}^N) + \frac{p_\eta^2}{2Q} + Nf k_B T_0 \eta
(15)

Where NfNf is the number of degrees of freedom.

Equations of Motion:

r˙i=pimi\dot{\mathbf{r}}_i = \frac{\mathbf{p}_i}{m_i}
(16)

p˙i=Fi−pηQpi\dot{\mathbf{p}}_i = \mathbf{F}_i - \frac{p_\eta}{Q}\mathbf{p}_i
(17)

η˙=pηQ\dot{\eta} = \frac{p_\eta}{Q}
(18)

p˙η=∑ipi2mi−NfkBT0\dot{p}_\eta = \sum_i \frac{p_i^2}{m_i} - Nf k_B T_0
(19)

Physical Interpretation:

Pros:

Cons:

Choosing QQ:

Q=NfkBT0τT2Q = Nf k_B T_0 \tau_T^2
(20)

where τT\tau_T is the desired coupling time (typically 0.1-1.0 ps).

6.3.3 Langevin Dynamics

Method: Add both friction and random force terms to mimic collisions with solvent:

mir¨i=Fi−γmir˙i+Ri(t)m_i \ddot{\mathbf{r}}_i = \mathbf{F}_i - \gamma m_i \dot{\mathbf{r}}_i + \mathbf{R}_i(t)
(21)

Where:

Pros:

Cons:

6.3.4 Summary: Choosing a Thermostat

ThermostatEnsembleUse Case
None (NVE)MicrocanonicalTesting energy conservation
Velocity RescalingApproximate NVTEquilibration
Nosé-HooverCanonical (NVT)Production runs, free energy
LangevinCanonical (NVT)Implicit solvent, large systems

6.4 Barostats (NPT)

To control Pressure, we treat the Volume (VV) of the simulation box as a dynamic variable.

Instantaneous Pressure: The virial theorem gives us:

P=NkBTV+13V⟨∑iri⋅Fi⟩P = \frac{Nk_B T}{V} + \frac{1}{3V}\left\langle \sum_i \mathbf{r}_i \cdot \mathbf{F}_i \right\rangle
(22)

The first term is the kinetic contribution, the second is the virial (force × distance).

6.4.1 Berendsen Barostat

Method: Scale the box size and all coordinates:

μ=[1−ΔtτPβT(P0−P(t))]1/3\mu = \left[1 - \frac{\Delta t}{\tau_P}\beta_T(P_0 - P(t))\right]^{1/3}
(23)

Where:

Update:

V(t+Δt)=μ3V(t)V(t + \Delta t) = \mu^3 V(t)
(24)

ri(t+Δt)=μri(t+Δt)\mathbf{r}_i(t + \Delta t) = \mu \mathbf{r}_i(t + \Delta t)
(25)

Pros/Cons: Same as Berendsen thermostat - good for equilibration, but doesn’t produce true NPT ensemble.

6.4.2 Parrinello-Rahman Barostat

Method: The box itself becomes a dynamical matrix h\mathbf{h} with its own equation of motion:

h¨=VW−1(Pinternal−P0)\ddot{\mathbf{h}} = V W^{-1}(\mathbf{P}_{internal} - \mathbf{P}_0)
(26)

Where:

Pros:

Cons:

6.5 Combined Thermostat + Barostat (NPT)

For NPT simulations, we typically use:

Practical Workflow:

  1. Energy minimization (eliminate bad contacts)

  2. NVT equilibration with Berendsen (heat up system)

  3. NPT equilibration with Berendsen (adjust density)

  4. NPT production with Nosé-Hoover/Parrinello-Rahman (data collection)

7. Statistical Mechanics Connection: The Ergodic Hypothesis

9.1 The Fundamental Question

If we simulate only 1000 atoms for 100 nanoseconds, how can we predict macroscopic properties of a mole (6.02×10236.02 \times 10^{23}) of atoms?

The solution is the Ergodic Hypothesis. It states that:

Over sufficiently long time periods, a single trajectory will visit all accessible microstates with the correct probability.

Therefore, the Time Average (from our simulation trajectory) equals the Ensemble Average (measuring all possible states simultaneously).

⟨A⟩time=lim⁡T→∞1T∫0TA(t)dt=⟨A⟩ensemble=∫A(rN,pN)ρ(rN,pN)drNdpN\langle A \rangle_{time} = \lim_{T \to \infty} \frac{1}{T}\int_0^T A(t) dt = \langle A \rangle_{ensemble} = \int A(\mathbf{r}^N, \mathbf{p}^N) \rho(\mathbf{r}^N, \mathbf{p}^N) d\mathbf{r}^N d\mathbf{p}^N
(27)

9.2 Practical Implications

What this means:

Caveats:

9.3 Statistical Inefficiency and Correlation

Not all data points are independent. If we save configurations every 1 fs, consecutive frames are highly correlated.

Autocorrelation Function:

C(t)=⟨A(t0)A(t0+t)⟩−⟨A⟩2⟨A2⟩−⟨A⟩2C(t) = \frac{\langle A(t_0) A(t_0 + t)\rangle - \langle A \rangle^2}{\langle A^2 \rangle - \langle A \rangle^2}
(28)

Correlation Time (τc\tau_c): The time for C(t)C(t) to decay to 1/e1/e of its initial value.

Statistical Inefficiency:

s=1+2∑t=1∞C(t)s = 1 + 2\sum_{t=1}^{\infty} C(t)
(29)

Effective Number of Independent Samples:

Neff=NtotalsN_{eff} = \frac{N_{total}}{s}
(30)

Practical Rule: Save configurations every 10τc10\tau_c to ensure independence.

8. Periodic Boundary Conditions

8.1 The Surface Problem

If we simulate a box of 1000 atoms with free boundaries, ~30% of atoms are at the surface. Surface atoms behave differently than bulk atoms.

Solution: Periodic Boundary Conditions (PBC)

The simulation box is replicated infinitely in all directions. When an atom leaves one side of the box, its periodic image enters from the opposite side.

8.2 Implementation

Minimum Image Convention: Each atom interacts with the nearest image of every other atom (including its own periodic images if the box is small).

Wrapping Coordinates:

def apply_pbc(positions, box_length):
    """Wrap positions into the central box"""
    return positions - box_length * np.floor(positions / box_length)

Computing Distances:

def minimum_image_distance(r1, r2, box_length):
    """Compute minimum image distance between two atoms"""
    dr = r1 - r2
    dr = dr - box_length * np.round(dr / box_length)
    return np.linalg.norm(dr)

8.3 Consequences

Advantages:

Disadvantages:

9. Summary and Workflow

9.1 Step-by-Step Protocol

1. System Preparation

2. Energy Minimization

3. Equilibration Phase

4. Production Phase

5. Analysis

9.2 Common Pitfalls and Best Practices

❌ Common Mistakes:

✅ Best Practices:

9.3 The Big Picture

When you run an MD simulation, you are orchestrating these components:

ComponentPhysical MeaningKey Parameters
Force FieldPotential Energy Surface U\mathcal{U}ϵ,σ\epsilon, \sigma (LJ); bond/angle constants
IntegratorTime Evolution (Verlet)Timestep Δt\Delta t (typically 1-2 fs)
ThermostatTemperature Control (NVT)Target T, coupling time τT\tau_T
BarostatPressure Control (NPT)Target P, coupling time τP\tau_P
BoundaryPeriodic ConditionsBox size (> 2×cutoff)
SamplingData CollectionOutput frequency

9.4 Topics to be Covered Later

The following advanced topics will be covered in detail in subsequent lectures:

10. Conclusion and Further Reading

10.1 What We’ve Learned

Molecular Dynamics is a powerful bridge between microscopic laws and macroscopic observations. By solving Newton’s equations for individual atoms, we can:

The key concepts covered in these notes are:

  1. Deterministic evolution in phase space via Hamilton’s equations

  2. Force fields that approximate the potential energy surface

  3. Symplectic integrators (Velocity Verlet) that conserve energy and phase space

  4. Thermostats/barostats that control temperature and pressure

  5. Ergodicity that connects time averages to ensemble averages

  6. Periodic boundary conditions that eliminate surface effects

These fundamental concepts form the foundation for all molecular dynamics simulations, from simple Lennard-Jones fluids to complex biomolecular systems.

Textbook:

10.3 The Future of MD

Current Frontiers:

The field continues to evolve rapidly, driven by increasing computational power and algorithmic innovations. The fundamental principles, however, remain those we’ve covered here.

End of Comprehensive Notes

Lecture Notes
Introduction to Statistical Mechanics
Lecture Notes
Practical Molecular Dynamics with LAMMPS