3 Molecular Dynamics

Molecular dynamics (MD) is a simulation technique in which we numerically integrate Newton’s equations of motion for a system of atoms. By collecting atomic trajectories, we can compute thermodynamic and dynamic properties via the ergodic hypothesis (Section 2.3.1 in Chapter 1).

3.1 Equations of Motion

For a system of $N$ atoms, Newton’s second law gives:

\[ m_i \ddot{\mathbf{r}}_i = \mathbf{F}_i = -\nabla_{\mathbf{r}_i} U(\mathbf{r}_1, \ldots, \mathbf{r}_N) \tag{3.1}\]

where $\mathbf{F}_i$ is the force on atom $i$ and $U$ is the potential energy surface. These $3N$ coupled second-order ODEs are what we integrate in MD.

The equations of motion conserve the total energy $\mathcal{H}$ (Eq. 2.2), so a perfect integrator would sample the NVE (microcanonical) ensemble.

3.2 Numerical Integration

3.2.1 The Velocity Verlet Algorithm

The most widely used integrator in MD is the velocity Verlet algorithm. Given positions $\mathbf{r}(t)$, velocities $\mathbf{v}(t)$, and forces $\mathbf{F}(t)$ at time $t$, the update to time $t + \Delta t$ is:

\[ \mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\,\Delta t + \frac{\mathbf{F}(t)}{2m}\,\Delta t^2 \tag{3.2}\]

\[ \mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{\mathbf{F}(t) + \mathbf{F}(t+\Delta t)}{2m}\,\Delta t \tag{3.3}\]

Note that Eq. 3.3 requires computing forces at the new positions before updating velocities — this is what distinguishes velocity Verlet from the simpler (less accurate) Euler method.

Choosing the timestep $\Delta t$

A good rule of thumb is $\Delta t \approx \tau_{\min}/20$, where $\tau_{\min}$ is the period of the fastest vibration in the system. For systems containing hydrogen, $\Delta t \sim 0.5$–$1$ fs. For heavy-atom-only systems (e.g., metals), $\Delta t \sim 2$–$5$ fs.

3.2.2 Conservation Properties

Velocity Verlet is symplectic — it conserves a shadow Hamiltonian close to the true $\mathcal{H}$, which prevents energy drift over long simulations. The total energy should fluctuate but not drift; a systematic drift indicates too large a timestep.

3.3 Temperature and Thermostats

The instantaneous kinetic temperature is defined via the equipartition theorem:

\[ T_{\text{inst}} = \frac{2 E_{\text{kin}}}{(3N - N_c)\, k_B} = \frac{1}{(3N - N_c)\,k_B} \sum_{i=1}^{N} m_i \mathbf{v}_i^2 \tag{3.4}\]

where $N_c$ is the number of constraints (e.g., 3 for fixed center-of-mass momentum).

To sample the NVT ensemble, we couple the system to a thermostat. The two most common thermostats are:

Table 3.1: Common thermostats and their properties.

Thermostat	Ensemble	Notes
Velocity rescaling	NVT (approx.)	Simple but non-rigorous
Nosé–Hoover	NVT	Rigorous, may have ergodicity issues
Nosé–Hoover chains	NVT	Most robust for most systems
Langevin	NVT	Stochastic; good for biomolecules

3.3.1 Nosé–Hoover Thermostat

The Nosé–Hoover equations of motion introduce a fictitious variable $\xi$ (the thermostat degree of freedom) with coupling parameter $Q$:

\[ \dot{\mathbf{r}}_i = \frac{\mathbf{p}_i}{m_i} \tag{3.5}\]

\[ \dot{\mathbf{p}}_i = \mathbf{F}_i - \xi\, \mathbf{p}_i \tag{3.6}\]

\[ \dot{\xi} = \frac{1}{Q}\left(\sum_i \frac{\mathbf{p}_i^2}{m_i} - (3N-N_c)\,k_BT\right) \tag{3.7}\]

When the kinetic energy exceeds the target $(3N-N_c)k_BT/2$, $\xi > 0$ and the equations of motion slow the particles down, and vice versa.

3.4 A Simple MD Simulation with ASE

Below is a minimal Python script using ASE to run an NVT MD simulation of bulk copper:

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.md.nvtberendsen import NVTBerendsen
from ase import units
from ase.io.trajectory import Trajectory

# Build the system
atoms = bulk('Cu', 'fcc', a=3.61) * (4, 4, 4)  # 4x4x4 supercell
atoms.calc = EMT()                               # Effective medium potential

# Set initial velocities at 300 K
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
MaxwellBoltzmannDistribution(atoms, temperature_K=300)

# Set up NVT MD (Berendsen thermostat)
dt = 5 * units.fs          # 5 fs timestep
md = NVTBerendsen(
    atoms,
    timestep=dt,
    temperature_K=300,
    taut=100 * units.fs,   # thermostat coupling time
    logfile='md.log',
    loginterval=10,
)

# Attach a trajectory writer
traj = Trajectory('cu_nvt.traj', 'w', atoms)
md.attach(traj.write, interval=10)

# Run 1000 steps
md.run(1000)
print("Simulation complete.")

3.5 Computing Properties from Trajectories

3.5.1 Radial Distribution Function

The radial distribution function (RDF) $g(r)$ quantifies the probability of finding an atom at distance $r$ from a reference atom, relative to an ideal gas:

\[ g(r) = \frac{V}{N^2} \left\langle \sum_{i \neq j} \delta(r - |\mathbf{r}_i - \mathbf{r}_j|) \right\rangle \tag{3.8}\]

3.5.2 Mean Squared Displacement and Diffusion

The diffusion coefficient $D$ is related to the mean squared displacement (MSD) via the Einstein relation:

\[ D = \lim_{t\to\infty} \frac{\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2 \rangle}{6t} \tag{3.9}\]

Cross-references in this chapter

Forces $\mathbf{F}_i$ in Eq. 3.1 come from the potential energy surface, discussed in Section 4.1 (Chapter 3).
Temperature Eq. 3.4 connects to the canonical ensemble (Section 2.2.1, Chapter 1).
The NVE ensemble sampled by bare velocity Verlet (Section 3.2.1) was defined in Table 2.1.

Exercises

Derive Eq. 3.2 and Eq. 3.3 from a Taylor expansion of $\mathbf{r}(t \pm \Delta t)$. What is the local truncation error in position?
Implement a simple 1D harmonic oscillator MD in Python using the velocity Verlet algorithm (Eq. 3.2, Eq. 3.3). Plot the total energy as a function of time for three different timesteps. At what timestep does the simulation become unstable?
Using Eq. 3.4, what is the kinetic temperature of a 64-atom copper system where the total kinetic energy is $0.5$ eV? ($k_B = 8.617 \times 10^{-5}$ eV/K)

--- title: "Molecular Dynamics" --- Molecular dynamics (MD) is a simulation technique in which we numerically integrate Newton's equations of motion for a system of atoms. By collecting atomic trajectories, we can compute thermodynamic and dynamic properties via the ergodic hypothesis (@sec-ergodic in Chapter 1). ## Equations of Motion {#sec-eom} For a system of $N$ atoms, Newton's second law gives: $$ m_i \ddot{\mathbf{r}}_i = \mathbf{F}_i = -\nabla_{\mathbf{r}_i} U(\mathbf{r}_1, \ldots, \mathbf{r}_N) $$ {#eq-newton} where $\mathbf{F}_i$ is the force on atom $i$ and $U$ is the potential energy surface. These $3N$ coupled second-order ODEs are what we integrate in MD. The equations of motion conserve the total energy $\mathcal{H}$ (@eq-hamiltonian), so a perfect integrator would sample the **NVE (microcanonical) ensemble**. ## Numerical Integration {#sec-integration} ### The Velocity Verlet Algorithm {#sec-verlet} The most widely used integrator in MD is the **velocity Verlet** algorithm. Given positions $\mathbf{r}(t)$, velocities $\mathbf{v}(t)$, and forces $\mathbf{F}(t)$ at time $t$, the update to time $t + \Delta t$ is: $$ \mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\,\Delta t + \frac{\mathbf{F}(t)}{2m}\,\Delta t^2 $$ {#eq-verlet-r} $$ \mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{\mathbf{F}(t) + \mathbf{F}(t+\Delta t)}{2m}\,\Delta t $$ {#eq-verlet-v} Note that @eq-verlet-v requires computing forces at the new positions before updating velocities — this is what distinguishes velocity Verlet from the simpler (less accurate) Euler method. ::: {.callout-tip} ## Choosing the timestep $\Delta t$ A good rule of thumb is $\Delta t \approx \tau_{\min}/20$, where $\tau_{\min}$ is the period of the fastest vibration in the system. For systems containing hydrogen, $\Delta t \sim 0.5$–$1$ fs. For heavy-atom-only systems (e.g., metals), $\Delta t \sim 2$–$5$ fs. ::: ### Conservation Properties Velocity Verlet is **symplectic** — it conserves a shadow Hamiltonian close to the true $\mathcal{H}$, which prevents energy drift over long simulations. The total energy should fluctuate but not drift; a systematic drift indicates too large a timestep. ## Temperature and Thermostats {#sec-thermostats} The instantaneous kinetic temperature is defined via the equipartition theorem: $$ T_{\text{inst}} = \frac{2 E_{\text{kin}}}{(3N - N_c)\, k_B} = \frac{1}{(3N - N_c)\,k_B} \sum_{i=1}^{N} m_i \mathbf{v}_i^2 $$ {#eq-temperature} where $N_c$ is the number of constraints (e.g., 3 for fixed center-of-mass momentum). To sample the NVT ensemble, we couple the system to a **thermostat**. The two most common thermostats are: | Thermostat | Ensemble | Notes | |------------|----------|-------| | Velocity rescaling | NVT (approx.) | Simple but non-rigorous | | Nosé–Hoover | NVT | Rigorous, may have ergodicity issues | | Nosé–Hoover chains | NVT | Most robust for most systems | | Langevin | NVT | Stochastic; good for biomolecules | : Common thermostats and their properties. {#tbl-thermostats} ### Nosé–Hoover Thermostat {#sec-nose-hoover} The Nosé–Hoover equations of motion introduce a fictitious variable $\xi$ (the thermostat degree of freedom) with coupling parameter $Q$: $$ \dot{\mathbf{r}}_i = \frac{\mathbf{p}_i}{m_i} $$ {#eq-nh-r} $$ \dot{\mathbf{p}}_i = \mathbf{F}_i - \xi\, \mathbf{p}_i $$ {#eq-nh-p} $$ \dot{\xi} = \frac{1}{Q}\left(\sum_i \frac{\mathbf{p}_i^2}{m_i} - (3N-N_c)\,k_BT\right) $$ {#eq-nh-xi} When the kinetic energy exceeds the target $(3N-N_c)k_BT/2$, $\xi > 0$ and the equations of motion slow the particles down, and vice versa. ## A Simple MD Simulation with ASE {#sec-ase-example} Below is a minimal Python script using ASE to run an NVT MD simulation of bulk copper: ```python from ase.build import bulk from ase.calculators.emt import EMT from ase.md.nvtberendsen import NVTBerendsen from ase import units from ase.io.trajectory import Trajectory # Build the system atoms = bulk('Cu', 'fcc', a=3.61) * (4, 4, 4) # 4x4x4 supercell atoms.calc = EMT() # Effective medium potential # Set initial velocities at 300 K from ase.md.velocitydistribution import MaxwellBoltzmannDistribution MaxwellBoltzmannDistribution(atoms, temperature_K=300) # Set up NVT MD (Berendsen thermostat) dt = 5 * units.fs # 5 fs timestep md = NVTBerendsen( atoms, timestep=dt, temperature_K=300, taut=100 * units.fs, # thermostat coupling time logfile='md.log', loginterval=10, ) # Attach a trajectory writer traj = Trajectory('cu_nvt.traj', 'w', atoms) md.attach(traj.write, interval=10) # Run 1000 steps md.run(1000) print("Simulation complete.") ``` ## Computing Properties from Trajectories {#sec-observables} ### Radial Distribution Function The radial distribution function (RDF) $g(r)$ quantifies the probability of finding an atom at distance $r$ from a reference atom, relative to an ideal gas: $$ g(r) = \frac{V}{N^2} \left\langle \sum_{i \neq j} \delta(r - |\mathbf{r}_i - \mathbf{r}_j|) \right\rangle $$ {#eq-rdf} ### Mean Squared Displacement and Diffusion The diffusion coefficient $D$ is related to the **mean squared displacement** (MSD) via the Einstein relation: $$ D = \lim_{t\to\infty} \frac{\langle |\mathbf{r}(t) - \mathbf{r}(0)|^2 \rangle}{6t} $$ {#eq-diffusion} ::: {.callout-note} ## Cross-references in this chapter - Forces $\mathbf{F}_i$ in @eq-newton come from the potential energy surface, discussed in @sec-forcefields (Chapter 3). - Temperature @eq-temperature connects to the canonical ensemble (@sec-canonical, Chapter 1). - The NVE ensemble sampled by bare velocity Verlet (@sec-verlet) was defined in @tbl-ensembles. ::: ## Exercises {.unnumbered} 1. Derive @eq-verlet-r and @eq-verlet-v from a Taylor expansion of $\mathbf{r}(t \pm \Delta t)$. What is the local truncation error in position? 2. Implement a simple 1D harmonic oscillator MD in Python using the velocity Verlet algorithm (@eq-verlet-r, @eq-verlet-v). Plot the total energy as a function of time for three different timesteps. At what timestep does the simulation become unstable? 3. Using @eq-temperature, what is the kinetic temperature of a 64-atom copper system where the total kinetic energy is $0.5$ eV? ($k_B = 8.617 \times 10^{-5}$ eV/K)