Module 2: The Physics Behind MD#
Welcome to Module 2! Here, we lay the theoretical groundwork for Molecular Dynamics simulations. We’ll explore the fundamental physics principles, how we model interactions, handle long-range forces, integrate motion, and apply boundary conditions.
Lesson 1: Classical Mechanics & Newton’s Equations#
At the core of MD simulations lies Classical Mechanics, specifically Newton’s Second Law of Motion. For each atom \(i\) in our system with mass \(m_i\) and position vector \(\mathbf{r}_i\), the law states:
where:
\(\mathbf{F}_i\) is the total force acting on atom \(i\) due to all other atoms in the system.
\(\mathbf{a}_i\) is the acceleration of atom \(i\).
Recall that acceleration is the second derivative of position with respect to time, \(\mathbf{a}_i = \frac{d^2 \mathbf{r}_i}{dt^2}\), and the first derivative of velocity, \(\mathbf{a}_i = \frac{d\mathbf{v}_i}{dt}\). So, we can rewrite Newton’s second law as a second-order ordinary differential equation:
The force \(\mathbf{F}_i\) depends on the positions of all atoms in the system (\(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N\)).
The Goal of MD: If we know the positions and velocities of all atoms at a certain time \(t\), and we have a way to calculate the force \(\mathbf{F}_i\) on each atom, we can use Newton’s equations to predict the positions and velocities at a slightly later time \(t + \Delta t\). By repeating this process over many small time steps \(\Delta t\), we simulate the trajectory of the system over time.
This process of predicting the future state based on the current state is called numerical integration. We’ll explore specific algorithms for this in Lesson 2.3.
Lesson 2: Potential Energy Surfaces & Force Fields#
Where do the forces (\(\mathbf{F}_i\)) come from? In classical mechanics, forces often arise from a potential energy function, \(V(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N)\), which depends on the positions of all particles in the system. This function defines the Potential Energy Surface (PES).
The force on atom \(i\) is the negative gradient of the potential energy with respect to the coordinates of atom \(i\):
This means if we can define the potential energy \(V\), we can calculate the forces needed for Newton’s equations.
Force Fields: In molecular simulations, the potential energy function \(V\) is called a Force Field (FF). It’s an empirical model, meaning it’s based on fitting mathematical functions to experimental data and/or quantum mechanical calculations. Force fields approximate the true quantum mechanical interactions with simpler classical functions, making calculations feasible for large systems.
A typical force field decomposes the total potential energy into several terms:
Bonded Terms: Describe interactions between atoms connected by covalent bonds (bond stretching, angle bending, dihedral torsion). These are typically short-ranged.
Non-Bonded Terms: Describe interactions between atoms that are not directly bonded (or are separated by several bonds).
van der Waals interactions (Lennard-Jones potential):
This term decays relatively quickly (approximately as \(r^{-6}\) in the attractive tail).
Electrostatic interactions (Coulomb potential):
This term decays more slowly (as \(r^{-1}\)).
Parameterization: The constants in these equations are the force field parameters, derived to reproduce experimental or QM data.
Lesson 3: Integrating the Equations of Motion#
Now that we know how to get forces from a potential (\(F = -\nabla V\)), how do we use them to update positions and velocities over time? We need a numerical integration algorithm to solve Newton’s equations (\(m\ddot{\mathbf{r}} = \mathbf{F}\)).
A very simple approach is the Euler method, but it is generally not accurate or stable enough for MD. A much more popular and suitable algorithm is the Velocity Verlet algorithm.
Velocity Verlet Algorithm#
Given positions \(\mathbf{r}\), velocities \(\mathbf{v}\), accelerations \(\mathbf{a}\), mass \(m\), and time step \(\Delta t\):
Half-step velocity update
Position update
Compute forces and new acceleration
Half-step velocity completion
Notes:
This scheme is time-reversible and conserves energy well for small \(\Delta t\).
Choose \(\Delta t\) to resolve the fastest motions (typically around 1-2 fs in atomistic MD).
def velocity_verlet_step(r_old, v_old, mass, dt, calculate_forces):
"""Performs one Velocity Verlet step with half-step velocity updates."""
# 0. Compute force & acceleration at current position
force_old = calculate_forces(r_old)
a_old = force_old / mass
# 1. Half-step velocity update
v_half = v_old + 0.5 * a_old * dt
# 2. Full position update using half-step velocity
r_new = r_old + v_half * dt
# 3. Compute force & acceleration at new position
force_new = calculate_forces(r_new)
a_new = force_new / mass
# 4. Complete velocity update with new acceleration
v_new = v_half + 0.5 * a_new * dt
return r_new, v_new
Lesson 4: Statistical Mechanics Snippets#
MD simulations generate microscopic information (positions, velocities). Statistical Mechanics connects this to macroscopic thermodynamic properties.
Thermodynamic Ensembles: Collections of microscopic states under macroscopic constraints.
Microcanonical (NVE): Constant Number (N), Volume (V), Energy (E). Natural for basic Verlet.
Canonical (NVT): Constant N, V, Temperature (T). Requires a thermostat.
Isothermal-Isobaric (NPT): Constant N, Pressure (P), T. Requires a thermostat and barostat.
Temperature from kinetic energy
Total kinetic energy:
Temperature estimate:
Here, \(N_{\mathrm{df}}\) is the number of degrees of freedom (for example, \(3N\) for \(N\) unconstrained particles in 3D), and \(k_B\) is the Boltzmann constant.
Lesson 5: Boundary Conditions#
Simulating finite systems introduces unrealistic surface effects. To simulate bulk materials, we use Periodic Boundary Conditions (PBC).
Imagine the simulation box is replicated infinitely in all directions.
Forces include interactions with particles in neighboring image boxes.
When a particle exits one face, it re-enters the opposite face.
(Visual Aid Idea: Include a diagram showing a central box and its periodic images in 2D)
+-------+-------+-------+
| | | |
| Box-1 | Box-0 | Box+1 |
| Image | Real | Image |
+-------+-------+-------+
| | | |
| Box-1 | Box-0 | Box+1 |
| Image | Real | Image |
+-------+-------+-------+
Lesson 6: Minimum Image Convention (MIC) & PBC Visualization#
Minimum Image Convention (MIC): When calculating the interaction between two particles \(i\) and \(j\) under PBC, we must consider the interaction between particle \(i\) and the closest periodic image of particle \(j\).
Assume a box of lengths \(L_x, L_y, L_z\).
Calculate the raw separation vector \(\Delta \mathbf{r} = \mathbf{r}_j - \mathbf{r}_i = (\Delta x, \Delta y, \Delta z)\).
Adjust each component using the nearest image distance:
\(\Delta x' = \Delta x - L_x\times{round}(\Delta x / L_x)\)
(Similarly for \(y\) and \(z\). The
round()function finds the nearest integer.)
The distance used for force/potential calculation is \(r_{ij} = |\Delta \mathbf{r}'|\).
This ensures we don’t interact with a distant image when a closer one is available and avoids double-counting interactions (especially important with cutoffs).
Visualizing PBC: It’s helpful to visualize how a particle “wraps around” the box. Imagine a 2D box:
If a particle’s \(x\) coordinate becomes greater than \(L_x\), its new position might be \(x' = x - L_x\).
If \(x\) becomes less than \(0\), its new position might be \(x' = x + L_x\). (Modulo arithmetic
x % Lcan also be used carefully).
We will implement functions for MIC and coordinate wrapping in Module 3.
End of Module 2. We have covered the core physics: Newton’s laws, potential energy functions (force fields), the issues of cutoffs and long-range forces (introducing PME), the Verlet algorithm for integration, the connection to thermodynamics (ensembles, temperature), and the concepts of Periodic Boundary Conditions (PBC) and the Minimum Image Convention (MIC).