MODULE 5

Molecular Dynamics Simulations

Introduction to Computational Chemistry (ICC)
Comprehensive Exam Notes

📋 Syllabus Coverage: Why Not Quantum Mechanics → Molecular Mechanics → Force Fields →
Bonded & Non-Bonded Terms → Equation of Motion → Integration Algorithms (Verlet) → MD
Algorithm → Simulation Cell → Periodic Boundary Conditions → Minimum Image Convention → Initial
Positions

1. Why Not Quantum Mechanics?

Ideally, to model complex molecular systems accurately, one would solve the full Schrödinger equation
for all subatomic particles.
(-ℏ²/2m) ∇²Ψ + U(x,y,z)Ψ(x,y,z) = EΨ(x,y,z)

Why QM Is Impractical for Large Systems

• Extremely difficult to program for complex multi-atom molecules
• Requires far more computing power than currently available for any system beyond a few atoms
• Scales exponentially with system size — computationally intractable
Key Insight: For most biochemical and chemical systems of interest (proteins, nucleic acids,
membranes), QM is accurate but computationally impossible. We use classical mechanics instead.

2. Classical / Molecular Mechanics (MM)

Instead of quantum mechanics, MD uses classical Newtonian mechanics to model molecular systems.

Core Principles of MM
• Calculates structure and energy of molecules based on nuclear motions only
• Electrons are NOT considered explicitly
• Based on the Born-Oppenheimer approximation of the Schrödinger equation:
◦ Nuclei are much heavier and move much slower than electrons
◦ Nuclear motions (vibrations, rotations) can be studied independently
◦ Electrons move fast enough to instantaneously adjust to any nuclear movement

The Ball-and-Spring Model

Molecular mechanics treats a molecule as a collection of weights connected by springs:
 Component Physical Representation Examples

Weights (balls) Atomic nuclei C, N, O, H atoms

Springs Chemical bonds C-C bond, C=O bond, N-H

bond

Three Pillars of Molecular Mechanics

• Atom Types — define each atom in a molecule (element + hybridization + environment)
• Parameters — numerical values for bond lengths, bond angles, force constants, etc.
• Equations — mathematical functions to calculate the energy of a molecule

3. MM: Atom Types

In MM, a given element may have several different atom types, depending on its chemical environment
and hybridization.

Why Multiple Atom Types?

• sp³-hybridized carbons have a tetrahedral bonding geometry → different type from aromatic C
• Aromatic carbons have a trigonal bonding geometry
• A C-C bond in an ethyl group differs from a C-C bond in a phenyl ring
• A C-C bond between a phenyl ring and an ethyl group differs from all other C-C bonds in
ethylbenzene
Each atom type carries its own mass, charge, van der Waals radius, and force constant parameters.
This fine-grained classification is what makes MM force fields accurate.

Typical Atom Type Properties

Property Description Example Value (Alanine N)

Charge Partial electrostatic charge on atom -0.47 e

Mass Atomic mass 14 (nitrogen)

Radius Van der Waals radius 1.7 Å

4. Force Field
4.1 Overview
Force Field: A complete set of atom types, parameters, and energy equations used to calculate the
energy and geometry of a molecule.
The total potential energy is split into bonded and non-bonded contributions:
V(R) = E_bonded + E_non-bonded
Force potentials are calculated independently and summed to give total energy. These equations
together define the Potential Energy Surface (PES) of the molecule.

4.2 Bonded Terms

Bonded terms describe interactions between atoms that are directly connected by chemical bonds.

(a) Bond Stretching

Models the energy required to stretch or compress a bond from its equilibrium length. Uses Hooke's
Law (harmonic approximation):
E_str = (1/2) k_str,ij (r_ij − r₀)²

Symbol Meaning

k_str,ij Force constant for stretching bond between atoms i and j (kcal/mol·Å²)

r_ij Current bond length between atoms i and j

r₀ Equilibrium (reference) bond length

Energy is zero at equilibrium (r_ij = r₀) and increases parabolically as the bond is stretched or
compressed.

(b) Bond Angle Bending

Models the energy cost of deforming a bond angle from its equilibrium value:
E_bend = (1/2) k_b,ijk (θ_ijk − θ₀)²

Symbol Meaning

k_b,ijk Force constant for bending angle i-j-k

θ_ijk Current angle between atoms i, j, k

θ₀ Equilibrium bond angle

(c) Dihedral (Torsional) Angle

Describes the energy barrier to rotation around a bond. Uses a cosine (periodic) function to capture
multiple energy minima:
E_d = k_d,ijkl [1 + cos(nΦ − Φ_s)]

Symbol Meaning

k_d,ijkl Torsional force constant

n Periodicity (number of minima per rotation)

Φ Current dihedral angle

Φ_s Phase shift (offset angle)

Unlike bonds and angles, dihedral energy is periodic — the function repeats as the bond rotates 360°.
Multiple minima correspond to eclipsed/staggered conformations.
(d) Improper Dihedral
Maintains planarity of flat molecular groups (e.g., peptide bond, aromatic rings):
E_id = (1/2) k_id,ijkl (Φ_ijkl − Φ₀)²
Used to keep atoms in a plane. Without this term, planar groups like benzene or the peptide
backbone would distort during simulation.

4.3 Non-Bonded Terms

Non-bonded terms describe interactions between atoms that are NOT directly bonded. They apply to
atom pairs separated by more than 3 bonds.

(a) Van der Waals Interaction — Lennard-Jones Potential

Captures short-range repulsion and long-range dispersion (London) attraction:
E_LJ = 4ε [(σ/r_ij)¹² − (σ/r_ij)⁶]

Symbol Meaning
ε (epsilon) Depth of the potential well (strength of interaction, in kcal/mol)

σ (sigma) Distance at which E_LJ = 0 (related to atomic size)

r_ij Distance between atoms i and j

(σ/r)¹² term Short-range Pauli repulsion (atoms resisting overlap)

(σ/r)⁶ term Long-range dispersion attraction (London forces)

The LJ potential has a minimum (most favorable interaction) at r = 2^(1/6) σ ≈ 1.12σ. At shorter
distances, repulsion dominates; at longer distances, energy → 0.

(b) Electrostatic Interaction — Coulomb's Law

Models charge-charge interactions between partial atomic charges:
F = q₁q₂ / (4πε₀ r²)

Symbol Meaning

q₁, q₂ Partial charges on atoms 1 and 2

ε₀ Permittivity of free space

r Distance between the charges

4.4 Complete Force Field Equation

The total force field energy is the sum of all terms:
E = Σ_bonds k_b(b₀-b)² + Σ_angles k_θ(θ₀-θ)² + Σ_planes k_ω(ω₀-ω)²
+ Σ_dihed k_φ(1-cos(nφ-δ)) + Σ_pairs (A/r¹² + B/r⁶ + Q_iQ_j/r_ij)

• Total energy = sum of all bond, angle, dihedral, improper, vdW, and electrostatic terms
 • Each force potential calculated independently, then summed
• These equations collectively define the Potential Energy Surface (PES)

5. Equation of Motion
5.1 Newton's Second Law
MD simulations use F = ma (Newton's second law) to determine how atoms move:
F = m·a → a_i = F_i / m_i
From the force field, we calculate the force on each atom as the negative gradient of potential energy:
F_ij = −∇U_ij

Step-by-Step Logic
Step Description

1. Calculate Force For each atom, sum forces from all other atoms: F_i = Σ_{j≠i} F_ij

2. Get Acceleration Use Newton's 2nd Law: a_i = F_i / m_i

3. Update Velocity Integrate: v = dr/dt

4. Update Position Integrate: r(t + Δt) from current position and velocity

Key Output: Integration of equations of motion yields a trajectory — a time-series of positions,

velocities, and accelerations for every atom. From the trajectory, thermodynamic properties can be
averaged.

Force Between Multiple Atoms

For a system of N atoms, each atom i experiences forces from all other atoms j:
F_i = Σ_{j≠i} F_ij (e.g., F₁ = F₁₂ + F₁₃)

5.2 Kinematic Analogy

The same equations apply as classical kinematics (just at the atomic scale):
Variable Formula Meaning

Velocity v = dr/dt Rate of change of position

Acceleration a = dv/dt Rate of change of velocity

Position update Δs = v₀t + ½at² Displacement under constant acceleration

6. Integration Algorithms
6.1 Why Numerical Integration?
The potential energy is a function of 3N atomic positions (for N atoms). Due to the complicated nature
of this function:
• There is NO analytical solution to the equations of motion
• Equations must be solved numerically using finite time steps
• Many algorithms have been developed, balancing accuracy vs. computational cost

6.2 The Verlet Algorithm (Original)

Derived from Taylor expansions of position forward and backward in time:
r(t + Δt) = r(t) + v(t)Δt + [f(t)/2m]Δt² + O(Δt⁴)
r(t − Δt) = r(t) − v(t)Δt + [f(t)/2m]Δt² − O(Δt⁴)
Adding both equations cancels the velocity term — giving the Verlet update formula:

r(t + Δt) ≈ 2r(t) − r(t − Δt) + [f(t)/m]Δt²

Velocity is recovered separately from positions at adjacent time steps:
v(t) = [r(t + Δt) − r(t − Δt)] / 2Δt

What the Algorithm Needs

Required Input Purpose

r(t) Current position

r(t − Δt) Previous position

f(t) Current force (from force field)

m Atomic mass

Properties of Verlet Algorithm

• Time-reversible (symmetric in time)
• Conserves energy well over long simulations
• Local error O(Δt⁴) — very accurate
• Drawback: velocity not directly available at the same time step as position

6.3 Other Common Integrators

Algorithm Key Feature

Velocity Verlet Positions and velocities computed at same time step; widely used

Leap-Frog Velocities computed at half-steps; good energy conservation

All MD programs (AMBER, CHARMM, NAMD) implement at least one of these integrators. Velocity
Verlet is the most commonly used in modern software.

7. The MD Algorithm
The complete MD algorithm proceeds as follows:
Step Action Details

1 Set Initial Conditions Assign positions (r) and velocities (v) for all atoms

2 Calculate Forces Compute F_i = −∇U for each atom using the force field

3 Update Positions Apply Verlet: r(t+Δt) = 2r(t) − r(t−Δt) + (f/m)Δt²

4 Update Velocities v(t) = [r(t+Δt) − r(t−Δt)] / 2Δt

5 Advance Time t → t + Δt
6 Repeat Continue steps 2-5 for hundreds of thousands of steps

7 Analysis Calculate thermodynamic properties from trajectory

The main output of an MD simulation is the TRAJECTORY — a record of all atomic positions (and
velocities) as a function of time. Average properties are computed from this trajectory.
Time step Δt ∼ 1–5 fs (femtoseconds, 10⁻¹⁵ s) — must be small enough to capture the fastest
atomic motions (bond vibrations).

8. Molecular Dynamics Simulations

8.1 What is MD?
MD simulations model the motion of atoms by solving classical equations of motion for a group of
particles.
MD = Force Field + Newton's Laws + Numerical Integration
MD is described as a 'computational microscope' — it reveals the molecular-level behavior of systems,
from electronic structures to long-timescale phase changes.

What MD Simulates
• Build a realistic atomistic model of the system
• Simulate system behavior over time under specific conditions (T, P, V, etc.)
• Track how atoms move, collide, and interact at femtosecond resolution

8.2 MD vs. QM Comparison

Aspect Quantum Mechanics (QM) Molecular Dynamics (MM/MD)

Electrons Explicitly treated Implicit (Born-Oppenheimer)

 Aspect Quantum Mechanics (QM) Molecular Dynamics (MM/MD)

Accuracy Higher Lower (approximation)

Speed Very slow Fast

System size Tens of atoms Millions of atoms

Timescale Femtoseconds Microseconds possible

9. The Simulation Cell

9.1 Setting Up the Box
Atoms are placed inside a 'porous' simulation super-cell (a box). The box size L is determined from the
experimental density:
L³ = M/D_exp = N·m / D_exp

Symbol Meaning

L Box side length

M Total mass of atoms in box

D_exp Experimental density

N Number of atoms

m Atomic mass per atom

Initial conditions: Positions are assigned from experiment (X-ray, NMR) or prior simulation. Initial
velocities are set to zero or drawn from a Maxwell-Boltzmann distribution.

9.2 Periodic Boundary Conditions (PBC)

The Problem PBC Solves
• A finite simulation box has surface atoms — atoms at the edges have fewer neighbors
• This creates artificial surface effects that distort bulk properties
• Solution: Replicate the simulation box infinitely in all directions

How PBC Works

• The simulation box is replicated to form an infinite periodic lattice
• Periodic images of each atom move with the same orientation and velocity
• When an atom leaves through one face, its image enters through the opposite face
• In 3D, the central box is surrounded by 26 immediate image cells

Image Cell Coordinates

x' = x + n₁L y' = y + n₂L z' = z + n₃L
where n₁, n₂, n₃ ∈ {-1, 0, 1} for the 26 nearest image cells.

Key Properties of PBC

Property Details

No surface atoms All atoms experience bulk-like environment

Linear momentum Conserved under PBC

Long-range interactions Require special treatment (Ewald summation) due to infinite re-
summation

Short-range interactions Handled using minimum image convention

9.3 Minimum Image Convention

Purpose
Prevents a particle from interacting with both the original particle j AND its image — which would
double-count interactions.

The Rule
• For any pair (i, j), only the CLOSEST image of j interacts with i
• Interactions beyond a cutoff distance Rc are ignored: U(Rc) ≈ 0
• Rc is chosen such that Rc ≤ L/2 (ensures each particle sees only one image)
Minimum Image Convention Rule: Particle i interacts with particle j OR one of j's images, but NEVER
both. The system size must satisfy Rc ≤ L/2.

Summary of PBC + Minimum Image

Concept Description

Simulation super-cell The primary box containing N atoms

Image cells Copies of the primary box tiling 3D space

PBC Atoms exiting one face re-enter from the opposite face

Minimum image convention Use closest image of each neighbor for force calculation

Cutoff Rc Maximum distance for non-bonded interactions; Rc ≤ L/2

10. Initial Positions

10.1 Sources of Initial Positions
Initial atomic coordinates (positions) for MD simulations can be obtained from:
Source Description Common Use

X-ray Crystallography High-resolution 3D structure from Proteins, nucleic acids

 Source Description Common Use

crystal diffraction

NMR Spectroscopy Solution-phase structure from Smaller proteins, flexible systems

nuclear magnetic resonance

Prior Simulation Output coordinates from a previous Continuation/restart of simulation

MD run

Analytical Calculation Computationally generated Small molecules, ideal structures

coordinates

10.2 The PDB Format

Atomic coordinates are stored in PDB (Protein Data Bank) format — a standard text-based format
where each line represents one atom:
ATOM 1 N NMETA 1 14.800 29.010 30.810 1.00 0.00

PDB Field Information

Record type ATOM or HETATM

Atom number Sequential atom index

Atom name Chemical identity + position (e.g., CA = alpha carbon)

Residue name Amino acid or molecule (e.g., ALA, NMETA)

Residue number Position in the chain

X, Y, Z coordinates Position in Angstroms (the highlighted columns in red)

Occupancy Fraction of time atom is in this position (usually 1.00)

B-factor Thermal fluctuation measure (0.00 for ideal structures)

The Protein Data Bank ([Link]/pdb/) is the primary repository for experimentally determined
macromolecular structures. It contains over 200,000 structures!

11. Quick Reference: Key Equations

Equation Formula Used For

Schrödinger (QM) (-ℏ²/2m)∇²Ψ + UΨ = EΨ Exact but impractical for large

systems

Bond Stretch E = ½k(r - r₀)² Bonded term: bond length deviation

Angle Bend E = ½k(θ - θ₀)² Bonded term: angle deviation

Dihedral E = k[1 + cos(nΦ - Φₛ)] Bonded term: torsional rotation

Improper Dihedral E = ½k(Φ - Φ₀)² Planarity constraint

 Equation Formula Used For

Lennard-Jones E = 4ε[(σ/r)¹² - (σ/r)⁶] Non-bonded: vdW interaction

Coulomb F = q₁q₂/(4πε₀r²) Non-bonded: electrostatics

Newton 2nd Law a_i = F_i / m_i Atomic acceleration from force

Verlet Position r(t+Δt) = 2r(t) - r(t-Δt) + (f/m)Δt² Advance atom positions

Verlet Velocity v(t) = [r(t+Δt) - r(t-Δt)] / 2Δt Calculate atom velocities

Box size L³ = N·m / D_exp Simulation cell dimensions

Image coordinates x' = x + n₁L, etc. Periodic boundary conditions

EXAM TIP: Know the physical meaning of every symbol in every equation. Be able to identify which
term is bonded vs. non-bonded, and why specific functional forms (harmonic, cosine, LJ) are used.