Molecular Dynamics II

Potentials

360.243 NSSC II — TU Wien

Ass. Prof. Esther Heid · Dr. Nico Unglert

This slide deck is interactive and best viewed in a browser.
In case you have a pdf version of it, I recommend viewing the slides at
https://nunglert.github.io/contents/presentations/nssc2_md

Use arrow keys to navigate · Press F for fullscreen

Outline

I. Quick MD recap

Numerical integration
Periodic boundary conditions
Thermodynamic ensembles
Trajectory analysis

II. Classical Potentials

Two-body — LJ, Buckingham, Morse
Solid state — SW, Tersoff
Molecular mechanics — OPLS-AA
Cutoffs and combination rules

III. Beyond classical

Quantum-mechanical methods
Machine-learned potentials
Descriptors & message passing

Companion exercise — soft-sphere LJ MD in Python

Integration algorithms

Truncating the Taylor expansion: Euler vs Velocity Verlet

Euler

\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \Delta t\,\mathbf{v}(t) + \tfrac{1}{2}\Delta t^2\,\mathbf{a}(t)$$ $$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \Delta t\,\mathbf{a}(t)

Simple, but not symplectic — total energy drifts.

Velocity Verlet

\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \Delta t\,\mathbf{v}(t) + \tfrac{1}{2}\Delta t^2\,\mathbf{a}(t)$$ $$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \Delta t\,\frac{\mathbf{a}(t) + \mathbf{a}(t+\Delta t)}{2}

Three-stage practice: r from $\mathbf{r}, \mathbf{v}, \mathbf{a}$ $\to$ a from new positions $\to$ v from old + new a.

Energy drift: Euler vs Velocity Verlet

Both are 2nd-order in $\Delta t$, but Velocity Verlet is symplectic, it conserves a slightly modified energy exactly, so the true $E$ only oscillates around its initial value rather than drifting.

Δt 0.15 steps 250 Velocity Verlet Euler exact

NVE on a 2D potential energy surface

Click the landscape to launch a deterministic trajectory — no thermostat, energy is conserved

Δt 0.020 steps 500 click anywhere on the PES to relaunch

PBC and the minimum-image partner

Hover an atom to highlight its closest periodic image of every other atom

Cutoff status

—

What you're seeing

The bold square is the central cell; the eight pale squares are its periodic copies. The selected (red) atom is connected by dashed lines to the nearest image of every other atom — what the simulation actually uses.

Hover any atom to switch the selection. Push $R_c$ past $L/2$ and the cutoff disc overlaps its own periodic image — that's the constraint $R_c < L/2$ made visible.

L 5.0 R_c 2.00 Show cutoff disc Show ghost cells MIC box

Thermodynamic ensembles

Different macroscopic constraints — different statistics

NVE (microcanonical)

$N$, $V$, $E$ fixed. Equal probability for all microstates at energy $E$.

\rho(\mathbf{r},\mathbf{p}) = \frac{\delta\!\left(H(\mathbf{r},\mathbf{p})-E\right)}{\Omega(N,V,E)}

\Omega = \int \delta\!\left(H - E\right) d\mathbf{r}\,d\mathbf{p}

Plain Newtonian dynamics — nothing extra needed.

NVT (canonical)

$N$, $V$, $T$ fixed. System exchanges energy with a heat bath at $\beta = 1/k_BT$.

\rho(\mathbf{r},\mathbf{p}) = \frac{e^{-\beta H(\mathbf{r},\mathbf{p})}}{Z(N,V,T)}

Z = \int e^{-\beta H(\mathbf{r},\mathbf{p})}\,d\mathbf{r}\,d\mathbf{p}

Requires a thermostat (Nosé–Hoover, Langevin, …).

NPT (isobaric–isothermal)

$N$, $P$, $T$ fixed. Volume fluctuates; system exchanges both energy and work.

\rho(\mathbf{r},\mathbf{p},V) = \frac{e^{-\beta(H + PV)}}{\Delta(N,P,T)}

\Delta = \int e^{-\beta(H + PV)}\,d\mathbf{r}\,d\mathbf{p}\,dV

Requires thermostat + barostat (Parrinello–Rahman, …).

NVT on the same PES — Andersen thermostat in action

Velocity resampling at rate $\nu$ lets the trajectory cross barriers and explore the full landscape

T 0.30 ν 1.0 steps 500 ν = 0 reduces to NVE

Live demo: 2D soft-sphere LJ — $g(r)$ builds as the system runs

N = 64 atoms, velocity Verlet, periodic boundaries, velocity-rescaling thermostat

Live readouts

E_k/N: —

T(t): — target:

g(r) samples: 0

Low T, high ρ: sharp peaks in g(r) — the system orders.
High T: g(r) flattens toward 1 — gas-like.
The first peak sits near $r = 2^{1/6}\sigma \approx 1.12$, the soft-sphere "contact" distance.

T 0.60 ρ 0.85

Ergodicity: a long trajectory samples the Boltzmann distribution

Analytical $\varrho \propto e^{-\beta U}$ on the left, NVT MD samples on the right

T 0.30 steps 3000 low T → concentrated near global min; high T → broad sampling

Part II — Potentials

Computing $E_{\text{pot}}(\mathbf{x})$

Three families, each a different point on the cost / accuracy curve:

Quantum mechanics — accurate, slow (DFT, post-HF)
Classical force fields — simple analytic forms fit to experiment, fast
Machine-learned potentials — flexible regressors fit to QM data

Constraints on any $E_{\text{pot}}$

Translation- and rotation-invariant
Permutation-invariant in identical particles
Differentiable — we need $\mathbf{f} = -\nabla E_{\text{pot}}$

The canonical pair potential: Lennard-Jones

V_{\text{LJ}}(r) = 4\varepsilon\!\left[\!\left(\tfrac{\sigma}{r}\right)^{\!12} \!-\! \left(\tfrac{\sigma}{r}\right)^{\!6}\right]

Two parameters fully determine the shape:

$\varepsilon$ — well depth (energy scale)
$\sigma$ — the distance at which $V=0$ (size scale)

Minimum at $r_{\min} = 2^{1/6}\sigma \approx 1.12\,\sigma$, depth $-\varepsilon$.

Force

F(r) = -\frac{dV_\mathrm{LJ}}{dr} = \frac{24\varepsilon}{r}\!\left[2\!\left(\tfrac{\sigma}{r}\right)^{\!12} \!-\! \left(\tfrac{\sigma}{r}\right)^{\!6}\right]

Soft-sphere variant

Truncate & shift at $r_{\min} = 2^{1/6}\sigma$ → purely repulsive. Used in the exercise notebook.

ε 1.00 σ 1.00 Soft-sphere

n-body decomposition

Systematic expansion into pair, triplet, … contributions

\begin{aligned} E_{\text{pot}}(\{\mathbf r\}) =\;\; &V_0 \\ + \sum_I &V_1(\mathbf r_I) \quad+ \sum_{I,J} V_2(\mathbf r_I, \mathbf r_J) \\ + \sum_{I,J,K} &V_3(\mathbf r_I, \mathbf r_J, \mathbf r_K) \quad+ \sum_{I,J,K,L} V_4(\ldots) + \cdots \end{aligned}

Rarely extended beyond $V_4$ — higher orders rarely buy accuracy worth their cost
Each $V_n$ typically depends only on rotation- and translation-invariant scalars (distances, angles, dihedrals)
The next slides walk through 2-body (bond), 3-body (angle), 4-body (dihedral) one by one

Two-body interaction: bond term

Expand any pair potential around its minimum → harmonic bond

V(r) \approx \sum_{k=0}^{n} \frac{V^{(k)}(r_{\min})}{k!}\,(r - r_{\min})^k \;=\; \underbrace{-\varepsilon}_{\text{const}} + \underbrace{\tfrac{1}{2}V''(r_{\min})(r-r_{\min})^2}_{\text{harmonic (order 2)}} + \cdots

Reading the curves

Solid (blue) — exact Lennard-Jones $V(r)$.
Dashed (orange) — Taylor expansion up to order $n$.
Green dot marks $r_{\min} = 2^{1/6}\sigma$ — the expansion point.

Since $V'(r_{\min}) = 0$, the series starts at order 2. Higher orders extend the accurate range outward from the minimum.

order n 2 true LJ Taylor

Three-body interaction: angle term

V(\theta) = K_\theta\,(\theta - \theta_0)^2

θ 109° K_θ 1.00 θ₀ 109°

Four-body interaction: Dihedral term

V(\phi) = \tfrac{1}{2}K_1\,(1 + \cos\phi) + \tfrac{1}{2}K_2\,(1 - \cos 2\phi) + \tfrac{1}{2}K_3\,(1 + \cos 3\phi)

φ 180° K₁0.4 K₂0.0 K₃1.2 total V(φ) n=1 n=2 n=3

Two-body potential flavors

LJ, Morse, Buckingham — repulsive core + attractive well, different mathematical form

\text{LJ:}\quad 4\varepsilon\!\left[\!\left(\tfrac{\sigma}{r}\right)^{\!12} \!-\! \left(\tfrac{\sigma}{r}\right)^{\!6}\right]

\text{Morse:}\quad D_e\!\left[1 - e^{-\alpha(r-R_e)}\right]^{\!2} \!-\! D_e

\text{Buck.:}\quad A\,e^{-Br} - \tfrac{C}{r^6}

Sliders control LJ ($\varepsilon$, $\sigma$); Morse and Buckingham are auto-tuned to match the same minimum position and well depth — the difference is shape only.

LJ: algebraic $r^{-12}$ repulsion — cheap, no physical justification for the exponent
Morse: exponential repulsion, finite at $r{=}0$; tunable width $\alpha$
Buckingham: $e^{-Br}$ repulsion (more physical) + $r^{-6}$ dispersion; can collapse at small $r$

ε 1.00 σ 1.00 Morse Buckingham Soft-sphere

Solid-state and bond-order potentials

When the angle between three atoms matters

Stillinger-Weber (1985)

$V_2(r)$ plus a three-body term penalising deviations from the tetrahedral angle:

\begin{aligned} V_3 =\, &\lambda \varepsilon\, (\cos\theta_{ijk} - \cos\theta_0)^2 \\ &\cdot\, e^{\gamma\sigma / (r_{ij} - a\sigma)}\, e^{\gamma\sigma / (r_{ik} - a\sigma)} \end{aligned}

Captures Si, GaN, and 2D crystals by enforcing local tetrahedral geometry.

Tersoff (1988) — bond order

Attractive part modulated by bond order $b_{ij}$ from the local environment:

V = f_c(r_{ij})\big[f_R(r_{ij}) + b_{ij}\, f_A(r_{ij})\big]

neighbours N 3

Molecular mechanics

Molecular systems (with a bent toward organic molecules in fluid phases)

Potential has bonded and non-bonded interactions
Atoms have vdW radii, partial charges and possibly polarizabilities
Molecules can be treated in an all-atom, united-atom (implicit hydrogen) or coarse-grained (multiple atom form a bead) fashion. Parts of the simulation can also be at the QM-level (QM–MM hybrid approach)

all-atom (AA)

united-atom (UA)

coarse-grained (CG)

Molecular mechanics — OPLS-AA

A typical biomolecular force field: bonded + non-bonded

V = V_{\text{bonds}} + V_{\text{angles}} + V_{\text{dihedrals}} + V_{\text{non-bonded}}

Bonds: $V_{\text{bond}} = \sum_{IJ} K_R [r_{IJ} - r_0]^2$
Angles: $V_{\text{angle}} = \sum_{IJK} K_\theta [\theta_{IJK} - \theta_0]^2$
Dihedrals: $V_{\text{dih}} = \sum_n \tfrac{1}{2} K_n [1 \pm \cos(n\phi - \phi_n)]$ — multiple cosine components per torsion
Non-bonded: LJ + partial-charge Coulomb, with a 1–4 fudge factor $f_{IJ} \in \{0, 0.5, 1\}$ by bond separation

Combination rules $A_{IJ} = \sqrt{A_I A_J},\;C_{IJ} = \sqrt{C_I C_J}$ keep the parameter count manageable.

Potential cutoffs

Three flavours

Truncate: drop V to 0 at $R_c$ — jump in energy and force.
Shift: subtract $V(R_c)$ so $V$ is continuous — force still jumps.
Switch: multiply by a smooth taper from $R_s$ to $R_c$ — both V and force smooth.

Cubic switch: $S(r) = 1 - 3t^2 + 2t^3$ with $t = (r-R_s)/(R_c-R_s)$; $S'(R_s) = S'(R_c) = 0$.

R_c1.55 R_s1.35 raw truncate shift switch

Part III — Beyond classical

Quantum mechanics, machine-learned potentials

Basic concepts

1900 Planck: Quantization of light

E = h\nu = \hbar\omega

energy · constant · frequency

1905 Einstein: Quantization due to photons

E = pc \;(\text{photons}), \quad \omega = ck

momentum · speed of light · wavenumber $k=2\pi/\lambda$

1913 Bohr: Quantized electron orbits

1924 De-Broglie: Standing wave model

E=\hbar\omega,\quad p=\hbar k \quad\text{(also for particles)}

1926 Schrödinger: Wave equation

i\hbar\,\frac{\partial\Psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2\Psi + V(\mathbf{x})\Psi

or $E\Psi = -\tfrac{\hbar^2}{2m}\nabla^2\Psi + V(\mathbf{x})\Psi$

1926 Born: Wave = probability amplitude

|\Psi|^2 \;\cong\; \text{probability of finding a particle at a given location}

Classical vs quantum mechanics

Classical mechanics

Classical Hamiltonian ($H: \Gamma \mapsto \mathbb{R} $)

H(\mathbf r,\mathbf p) = \frac{|\mathbf p|^2}{2m} + V(\mathbf r)

A function eating a point in phase space $(r, p) \in \Gamma$

Equations of motion (Hamilton)

\dot{\mathbf r} = \frac{\partial H}{\partial \mathbf p},\qquad \dot{\mathbf p} = -\frac{\partial H}{\partial \mathbf r}

$H(\mathbf r,\mathbf p)$ is the conserved energy.

Quantum mechanics

Quantum Hamiltonian ($ \hat{H}: \mathcal{H} \mapsto \mathcal{H} $)

\hat H = -\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf r)

now an operator eating a wavefunction $\psi \in \mathcal{H}$, where $\mathcal{H}$ is a Hilbert space

Equation of motion (Schrödinger)

i\hbar\,\frac{\partial \psi}{\partial t} = \hat H\,\psi

Stationary states: $\hat H\phi = E\phi$

Simplest case: particle in a box

A single quantum particle confined to $[0, L]$ — standing-wave eigenstates

Boundary conditions

Infinite walls at $x = 0$ and $x = L$: $\phi(0) = \phi(L) = 0$. Inside, $V = 0$ so the Schrödinger equation reduces to $\phi'' = -\frac{2mE}{\hbar^2}\phi$.

Eigenstates

\phi_n(x) = \sqrt{\tfrac{2}{L}}\,\sin\!\left(\tfrac{n\pi x}{L}\right), \qquad n = 1, 2, 3, \ldots

Energies

E_n = \frac{n^2 \pi^2 \hbar^2}{2 m L^2}

Energy is quantised: discrete levels labelled by $n$, spacing growing like $n^2$.

n 1 show ψ(x) show |ψ|² show energy levels

What the Schrödinger equation gives us

Hydrogenic orbitals — xz cross-section of $\psi(\mathbf{r})$

For the hydrogen atom we can formulate the Hamiltonian in atomic units

\hat{H} = -\tfrac{1}{2}\nabla^2 - \tfrac{1}{r} \quad \hat{H} \psi = E \psi

Each orbital factorises as $\psi_{n\ell m}(r,\theta,\phi) = R_{n\ell}(r)\,Y_\ell^m(\theta,\phi)$. The 2D slice through $y = 0$ makes radial nodes ($R = 0$) and angular nodes ($Y = 0$) visible as white separators between red (+) and blue (−) lobes.

orbital zoom 1.0× show |ψ|²

Many-body electronic structure

Fix the nuclei (Born–Oppenheimer), then solve for $N_e$ interacting electrons

Electronic Hamiltonian

With nuclei fixed at $\{\mathbf R_I\}$, the Born–Oppenheimer Hamiltonian for $N_e$ electrons is

\hat H = -\sum_i \frac{\hbar^2}{2 m_e}\nabla_i^2 - \sum_{I,i} \frac{Z_I e^2}{|\mathbf r_i - \mathbf R_I|} + \sum_{i \lt j} \frac{e^2}{|\mathbf r_i - \mathbf r_j|}

Eigenstate problem:

\hat H\,\Psi(\mathbf r_1,\ldots,\mathbf r_{N_e}) = \textcolor{red}{E} \,\Psi(\mathbf r_1,\ldots,\mathbf r_{N_e})

Connection to Classical stat. mech.:

\textcolor{red}{E} \equiv E(\{\mathbf R_I\}) = E_\mathrm{pot}(\{ \mathbf x \})

Now using ${\mathbf x}$ for the nuclear positions, like on earlier slides

Molecules

HF / DFT

Mean field. $\mathcal{O}(N^3)$.

MP2/CCSD

Perturbative. $\mathcal{O}(N^5)$–$N^6$.

CCSD(T)

“Gold std.” $\mathcal{O}(N^7)$.

FCI

Exact. Exponential cost.

Solids

DFT

Workhorse.

QP self-energy. Band gaps.

RPA

Correlation via response.

DMFT

Strong correlations.

Machine-learned potentials

Two key assumptions underlie all ML force fields

1. Local additivity

E_{\text{pot}} \approx \sum_i E_i\!\left(\mathcal{N}_i(R_c)\right)

Select an atom — neighbours within $R_c$ contribute to its energy.

2. Invariant descriptor

E_i = E_\theta\!\left[p(\mathbf{x}_i)\right]

$p(\mathbf{x})$ encodes the local environment, invariant to translation & rotation.

R_c 2.5

Anatomy of a neural-network potential

Atomic decomposition: one shared network $E_\theta$ per atom, summed to $E_{\text{total}}$

Local energy: $E_\theta(p(\mathbf{x}_i))$

\mathbf{h}^{(l)} = f\!\left(\mathbf{W}^{(l)}\mathbf{h}^{(l-1)}+\mathbf{b}^{(l)}\right)

descriptor p 0.50

$E_{\text{total}}(\mathbf{x}) = \sum_i E_\theta\!\left(p(\mathbf{x}_i)\right)$ — one network $E_\theta$ applied to every atom's local descriptor $p(\mathbf{x}_i)$. Forces: $\mathbf{F}_i = -\nabla_{\mathbf{x}_i} E_{\text{total}}$.

Descriptor generation: Message passing (Graph-NNs)

Update rule

m_{ij}^{(t)} = \phi_m\!\left(h_i^{(t)},\, h_j^{(t)}\right)

h_i^{(t+1)} = \phi_u\!\Big(h_i^{(t)},\; \bigoplus_{j \in \mathcal{N}(i)} m_{ij}^{(t)}\Big)

$h_i^{(t)}$: atom state
$m_{ij}^{(t)}$: message $j{\to}i$
$\phi_m,\phi_u$: message/update function (learnable)
$\mathcal{N}(i)$: neighbours

Now the descriptor generation itself is trainable and hence depends on a set of parameters $\theta$, i.e. $p_\theta(\mathbf x_i)$

round t = 0 reached: — Click any atom to switch source

Do you want to learn more about…

Classical potentials
Advances over classical potentials: Polarizability, lone-pairs, reactive force-fields, cross terms
How to treat electrostatics with periodic boundary conditions
Quantum chemistry, density functional theory
State-of-the-art machine-learning potentials
AI in chemistry

165.144 VU Simulation of condensed matter: Lecture and exercises

165.166 VU Introduction to Theoretical Chemistry: Lecture and exercises

165.164 PR Selected Topics in Theoretical Chemistry Research project

165.171 VU AI in Chemistry: Lecture and exercises