NE 451: Simulation Methods

Conrard Tetsassi Feugmo

Estimated study time: 55 minutes

Table of contents

Sources and References

Primary references — D. Frenkel & B. Smit, Understanding Molecular Simulation, 2nd ed., Academic Press, 2002; D.S. Sholl & J.A. Steckel, Density Functional Theory: A Practical Introduction, Wiley, 2009; D.P. Landau & K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, 4th ed., Cambridge, 2014; M.P. Allen & D.J. Tildesley, Computer Simulation of Liquids, 2nd ed., Oxford, 2017.

Online resources — LAMMPS documentation (lammps.org); Quantum ESPRESSO documentation (quantum-espresso.org); J.M. Thijssen, Computational Physics, 2nd ed., Cambridge, 2007; NIST Interatomic Potentials Repository (ctcms.nist.gov/potentials); J. Behler, “Atom-centered symmetry functions for constructing high-dimensional neural network potentials,” J. Chem. Phys. 134, 074106 (2011); K. Schütt et al., “SchNet: A deep learning architecture for molecules and materials,” J. Chem. Phys. 148, 241722 (2018); S. Batzner et al., “E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials,” Nat. Commun. 13, 2453 (2022).

Chapter 1: Foundations — Classical Statistical Mechanics

1.1 Why Simulation?

Materials science and nanotechnology operate at length and time scales that span many orders of magnitude. At the nanometre scale, individual atoms and electrons determine bulk properties: the band gap of a semiconductor, the fracture toughness of a ceramic, the selectivity of a catalyst. Analytical solutions to the governing equations — Newton’s laws for nuclei, Schrödinger’s equation for electrons — are unavailable for systems of more than a handful of interacting particles. Computational simulation bridges the gap between first principles and experiment by constructing a tractable mathematical model of a physical system and evolving it either deterministically or stochastically.

Four methodologies dominate nanotechnology simulation. Molecular dynamics (MD) integrates Newton’s equations of motion for a classical many-body system, yielding dynamical information: trajectories, diffusion coefficients, viscosities, elastic constants. Monte Carlo (MC) samples configuration space according to the Boltzmann distribution without following physical trajectories, granting access to equilibrium thermodynamic averages efficiently. Density functional theory (DFT) solves the quantum-mechanical electronic structure problem within a mean-field approximation, providing ground-state energies, forces, electronic structure, and optical properties. Machine learning (ML) methods learn interatomic interactions from DFT data sets, delivering DFT-level accuracy at a fraction of the computational cost, enabling access to length and time scales inaccessible by ab initio means.

Understanding when to use each method is as important as knowing how to use it.

Method selection guide. Use DFT when quantum electronic effects govern the property of interest — reaction barriers, band gaps, charge transfer. Use classical MD when dynamics at 10–100 nm and nanosecond time scales are needed and a validated force field exists. Use MC when equilibrium properties suffice and dynamic trajectories are unnecessary — adsorption isotherms, phase diagrams, defect concentrations. Use kinetic Monte Carlo when rare events (surface diffusion hops, nucleation events) control the evolution and a rate catalogue is available. Use ML potentials when DFT accuracy is needed at MD time and length scales.

1.2 The N-Body Problem in Classical Statistical Mechanics

Consider a system of \( N \) classical particles, each with position \( \mathbf{r}_i \in \mathbb{R}^3 \) and momentum \( \mathbf{p}_i \in \mathbb{R}^3 \). The microstate is the point \( \Gamma = (\mathbf{r}^N, \mathbf{p}^N) \) in the \( 6N \)-dimensional phase space. The total energy is the Hamiltonian

\[ \mathcal{H}(\Gamma) = \sum_{i=1}^{N} \frac{\mathbf{p}_i^2}{2m_i} + U(\mathbf{r}^N), \]

where \( U(\mathbf{r}^N) \) is the potential energy surface encoding all interparticle interactions.

Hamilton’s equations of motion,

\[ \dot{\mathbf{r}}_i = \frac{\partial \mathcal{H}}{\partial \mathbf{p}_i} = \frac{\mathbf{p}_i}{m_i}, \qquad \dot{\mathbf{p}}_i = -\frac{\partial \mathcal{H}}{\partial \mathbf{r}_i} = -\nabla_{\mathbf{r}_i} U, \]

are equivalent to Newton’s second law. For a system of \( N = 10^3 \) argon atoms this represents 6000 coupled first-order ODEs. For \( N = 10^6 \) it becomes 6 million. The system is in principle integrable only for \( N \leq 2 \); for larger systems we resort to numerical integration.

1.3 Ensembles and the Ergodic Hypothesis

A macroscopic thermodynamic state is not a single microstate but a statistical ensemble — a distribution over all microstates consistent with macroscopic constraints. Three ensembles appear throughout this course.

The microcanonical ensemble (NVE) fixes particle number \( N \), volume \( V \), and total energy \( E \). Each accessible microstate has equal probability \( 1/\Omega \), where \( \Omega(N,V,E) \) is the density of states. Entropy is \( S = k_B \ln \Omega \).

The canonical ensemble (NVT) fixes \( N \), \( V \), and temperature \( T \). The probability of microstate \( \Gamma \) is

\[ P(\Gamma) = \frac{1}{Z} e^{-\beta \mathcal{H}(\Gamma)}, \qquad \beta = \frac{1}{k_B T}, \]

where the canonical partition function is

\[ Z(N,V,T) = \frac{1}{N! h^{3N}} \int d\Gamma \, e^{-\beta \mathcal{H}(\Gamma)}. \]

The Helmholtz free energy is \( F = -k_B T \ln Z \).

The grand canonical ensemble (\( \mu \)VT) allows particle exchange with a reservoir at chemical potential \( \mu \). The grand partition function is

\[ \Xi(\mu, V, T) = \sum_{N=0}^{\infty} \frac{e^{\beta \mu N}}{N! h^{3N}} \int d\Gamma_N \, e^{-\beta \mathcal{H}_N}. \]

Ergodic hypothesis. A dynamical system is ergodic if, for almost every initial condition, the time average of any integrable observable \( A \) equals the ensemble average:

\[ \langle A \rangle_{\text{time}} = \lim_{T\to\infty} \frac{1}{T} \int_0^T A(\Gamma(t))\, dt = \langle A \rangle_{\text{ensemble}} = \int d\Gamma\, A(\Gamma) P(\Gamma). \]

The ergodic hypothesis underpins all of molecular dynamics: we simulate a single trajectory and compute time averages, then interpret them as thermodynamic ensemble averages.

The hypothesis fails for systems trapped in metastable states (glasses, proteins with high barriers), motivating enhanced-sampling techniques. In Monte Carlo, ergodicity requires that the Markov chain can reach any configuration from any other in a finite number of steps.

1.4 Thermodynamic Observables as Phase-Space Averages

All measurable thermodynamic quantities can be expressed as ensemble averages. The internal energy is simply \( U = \langle \mathcal{H} \rangle \). Temperature in the canonical ensemble is related to the kinetic energy by the equipartition theorem:

\[ \left\langle \frac{\mathbf{p}_i^2}{2m_i} \right\rangle = \frac{3}{2} k_B T \implies T = \frac{2}{3 N k_B} \left\langle \sum_i \frac{\mathbf{p}_i^2}{2m_i} \right\rangle. \]

Pressure is given by the virial theorem, derived in Chapter 4. The heat capacity at constant volume is

\[ C_V = \frac{\partial U}{\partial T}\bigg|_V = \frac{\langle \mathcal{H}^2 \rangle - \langle \mathcal{H} \rangle^2}{k_B T^2}, \]

showing that thermodynamic response functions are related to fluctuations — a profound consequence of statistical mechanics directly accessible in simulation.

Chapter 2: Molecular Dynamics — Theory and Algorithms

2.1 The Classical MD Algorithm

The core idea of molecular dynamics is to follow the time evolution of a classical \( N \)-body system by numerically integrating Newton’s equations of motion. The algorithm in pseudocode is:

Initialise positions \( \mathbf{r}_i \) and velocities \( \mathbf{v}_i \).
Compute forces \( \mathbf{F}_i = -\nabla_{\mathbf{r}_i} U(\mathbf{r}^N) \).
Integrate equations of motion over time step \( \Delta t \) to obtain new positions and velocities.
Update forces.
Accumulate observables. Repeat from step 3.

The central challenge is step 2: for a pair potential \( U = \sum_{i

2.2 Time Integration Algorithms

2.2.1 The Verlet Algorithm

The simplest and most widely used integrator is derived from two Taylor expansions of the position about time \( t \):

\[ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \mathbf{v}_i(t)\Delta t + \frac{1}{2}\mathbf{a}_i(t)(\Delta t)^2 + \frac{1}{6}\dot{\mathbf{a}}_i(t)(\Delta t)^3 + O((\Delta t)^4), \]\[ \mathbf{r}_i(t - \Delta t) = \mathbf{r}_i(t) - \mathbf{v}_i(t)\Delta t + \frac{1}{2}\mathbf{a}_i(t)(\Delta t)^2 - \frac{1}{6}\dot{\mathbf{a}}_i(t)(\Delta t)^3 + O((\Delta t)^4). \]

Adding these two equations:

\[ \mathbf{r}_i(t + \Delta t) = 2\mathbf{r}_i(t) - \mathbf{r}_i(t - \Delta t) + \mathbf{a}_i(t)(\Delta t)^2 + O((\Delta t)^4). \]

This is the Verlet position update. The position error per step is \( O((\Delta t)^4) \), making Verlet a fourth-order accurate integrator in position. Velocities can be estimated from

\[ \mathbf{v}_i(t) = \frac{\mathbf{r}_i(t+\Delta t) - \mathbf{r}_i(t-\Delta t)}{2\Delta t} + O((\Delta t)^2), \]

but these are only second-order accurate. A more numerically convenient form is the velocity-Verlet algorithm.

2.2.2 Velocity-Verlet Algorithm

Derivation of velocity-Verlet from Taylor expansion.

Taylor-expanding position and velocity to second order:

\[ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \mathbf{v}_i(t)\Delta t + \frac{1}{2}\mathbf{a}_i(t)(\Delta t)^2 + O((\Delta t)^3), \tag{1} \]\[ \mathbf{v}_i(t + \Delta t) = \mathbf{v}_i(t) + \mathbf{a}_i(t)\Delta t + \frac{1}{2}\dot{\mathbf{a}}_i(t)(\Delta t)^2 + O((\Delta t)^3). \tag{2} \]

The acceleration \( \mathbf{a}_i(t+\Delta t) = \mathbf{F}_i(\mathbf{r}^N(t+\Delta t))/m_i \) is computed after the position update. Approximating \( \dot{\mathbf{a}}_i(t) \approx [\mathbf{a}_i(t+\Delta t) - \mathbf{a}_i(t)]/\Delta t \) and substituting into (2):

\[ \mathbf{v}_i(t+\Delta t) = \mathbf{v}_i(t) + \frac{1}{2}\left[\mathbf{a}_i(t) + \mathbf{a}_i(t+\Delta t)\right]\Delta t + O((\Delta t)^3). \]

This yields the velocity-Verlet update:

Step 1. Advance positions: \( \mathbf{r}_i(t+\Delta t) = \mathbf{r}_i(t) + \mathbf{v}_i(t)\Delta t + \tfrac{1}{2}\mathbf{a}_i(t)(\Delta t)^2 \)

Step 2. Compute \( \mathbf{a}_i(t+\Delta t) \) from forces at new positions.

Step 3. Advance velocities: \( \mathbf{v}_i(t+\Delta t) = \mathbf{v}_i(t) + \tfrac{1}{2}\left[\mathbf{a}_i(t) + \mathbf{a}_i(t+\Delta t)\right]\Delta t \)

Velocity-Verlet is symplectic (it preserves the phase-space volume element) and time-reversible. It is algebraically equivalent to the original Verlet algorithm but self-starting and numerically superior, making it the standard in modern MD codes.

2.2.3 Higher-Order Integrators and Time Step Selection

Gear predictor-corrector schemes use information from several previous time steps to achieve fifth- or sixth-order accuracy in position. They are not symplectic and therefore do not conserve energy exactly over long runs, limiting their utility for NVE simulations. They excel in specialised situations such as stiff constraints.

Choosing the time step. The time step \( \Delta t \) must be smaller than the fastest dynamical mode in the system. For typical materials:

Liquids (LJ argon near melting): \( \Delta t \approx 1\text{–}2 \) fs.
Solid metals (EAM): \( \Delta t \approx 1\text{–}5 \) fs.
Covalent systems (Tersoff Si): \( \Delta t \approx 0.5\text{–}1 \) fs.
Systems with O-H or C-H bonds (classical): \( \Delta t \approx 0.5 \) fs or constrain H positions.

A safe criterion is that the total energy in an NVE run should drift by less than \( 10^{-4} k_B T \) per atom per nanosecond. Energy conservation is the primary diagnostic for integrator quality.

2.3 Boundary Conditions

2.3.1 Periodic Boundary Conditions

Periodic boundary conditions (PBC). To simulate a bulk material without surface effects using a small simulation box, we tile space with identical copies of the simulation cell. Particle \( i \) in the primary cell at position \( \mathbf{r}_i \) has an infinite set of images at \( \mathbf{r}_i + n_1 \mathbf{a}_1 + n_2 \mathbf{a}_2 + n_3 \mathbf{a}_3 \), where \( \mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_3 \) are the cell vectors and \( n_1, n_2, n_3 \in \mathbb{Z} \). A particle exiting through one face re-enters through the opposite face.

2.3.2 Minimum Image Convention

When computing the pair interaction between particles \( i \) and \( j \), we use the minimum image convention: the distance is taken to the nearest image of \( j \) relative to \( i \). For an orthorhombic box with side lengths \( L_x, L_y, L_z \):

\[ \Delta x_{ij} \to \Delta x_{ij} - L_x \cdot \text{round}\!\left(\frac{\Delta x_{ij}}{L_x}\right), \]

and similarly for \( y \) and \( z \). The cutoff radius \( r_c \) must satisfy \( r_c < L/2 \) in each dimension to ensure each particle interacts with at most one image of every other particle.

2.3.3 Slab Geometry and Non-Periodic Directions

For surfaces, interfaces, or thin films, PBC are applied in two dimensions (parallel to the surface) while the third dimension is non-periodic. A vacuum gap of at least 10–15 Å above the surface prevents spurious interactions across the simulation boundary.

Chapter 3: Interatomic Potentials and Force Fields

3.1 The Role of the Potential Energy Surface

The accuracy and transferability of a classical MD simulation is limited entirely by the quality of the interatomic potential. A potential energy surface (PES) must reproduce the energetics, forces, and elastic properties of the material across the configurations sampled during the simulation.

Force field. A force field is a mathematical model for the potential energy \( U(\mathbf{r}^N) \) of a classical many-body system, consisting of a functional form and a set of parameters fit to experimental data or ab initio calculations. Force fields range from simple pair potentials to many-body reactive forms capable of describing bond breaking and formation.

3.2 Lennard-Jones Potential

The Lennard-Jones (LJ) potential is the simplest realistic model for van der Waals interactions between noble-gas atoms:

Lennard-Jones potential. The pair potential between atoms \( i \) and \( j \) separated by distance \( r \) is

\[ u_{\text{LJ}}(r) = 4\varepsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right], \]

where \( \varepsilon \) is the depth of the potential well and \( \sigma \) is the finite distance at which \( u_{\text{LJ}} = 0 \). The minimum of the well occurs at \( r_{\min} = 2^{1/6}\sigma \) with depth \( -\varepsilon \). The \( r^{-6} \) term describes London dispersion (attractive); the \( r^{-12} \) term is a computational convenience for the steeply repulsive Pauli exclusion interaction.

For argon: \( \varepsilon / k_B = 119.8 \) K, \( \sigma = 3.405 \) Å. The LJ potential is computationally trivial and physically transparent; it is the workhorse of method development and benchmarking.

Example 3.1: LJ argon simulation setup in LAMMPS.

A minimal LAMMPS input script for liquid argon in reduced units:

units        lj
atom_style   atomic
lattice      fcc 1.0
region       box block 0 5 0 5 0 5
create_box   1 box
create_atoms 1 box
mass         1 1.0
pair_style   lj/cut 2.5
pair_coeff   1 1 1.0 1.0 2.5
velocity     all create 1.0 4928459
fix          1 all nvt temp 1.0 1.0 0.1
timestep     0.005
thermo       100
run          10000

Here all quantities are in reduced LJ units: energy in \( \varepsilon \), length in \( \sigma \), time in \( \sigma\sqrt{m/\varepsilon} \). The cutoff \( r_c = 2.5\sigma \) is standard for LJ; at this distance \( u_{\text{LJ}} \approx -0.016\varepsilon \). A tail correction to energy and pressure accounts for the truncated long-range contribution.

3.3 Embedded Atom Method

The LJ potential treats every atom as electronically isolated. For metals, the cohesive energy depends on the local electron density — a collective effect that pair potentials cannot capture. The Embedded Atom Method (EAM) of Daw and Baskes (1984) accounts for this:

\[ U = \sum_i F_i\!\left(\bar{\rho}_i\right) + \frac{1}{2}\sum_{i \neq j} \phi_{ij}(r_{ij}), \]

where \( \bar{\rho}_i = \sum_{j \neq i} \rho_j(r_{ij}) \) is the host electron density at atom \( i \) from all neighbours, \( F_i(\bar{\rho}_i) \) is the embedding energy, and \( \phi_{ij}(r_{ij}) \) is a pair repulsion. EAM potentials successfully reproduce elastic constants, surface energies, stacking fault energies, and vacancy formation energies for FCC and BCC metals. The NIST repository provides extensively validated EAM files for Fe, Cu, Al, Ni, and their alloys.

3.4 Tersoff Potential for Covalent Systems

Silicon, carbon, and other covalent materials form directional bonds whose energy depends on the local coordination environment. The Tersoff potential captures this through a bond-order formalism:

\[ U = \sum_{i} \sum_{j > i} f_c(r_{ij}) \left[ A e^{-\lambda_1 r_{ij}} - B_{ij} e^{-\lambda_2 r_{ij}} \right], \]

where \( B_{ij} \) is a many-body bond-order parameter that decreases as the coordination number increases, penalising over-coordination. The cutoff function \( f_c(r) \) smoothly transitions from 1 to 0. Tersoff parameters are available for Si, Ge, C, and SiC. The potential reproduces the diamond cubic structure, elastic constants, and phonon dispersion of silicon with qualitative accuracy.

3.5 Reactive Force Fields — ReaxFF

Both EAM and Tersoff are non-reactive: they cannot describe bond formation or breaking in response to changes in chemistry. ReaxFF (van Duin et al., 2001) uses a bond-order formalism where bond orders are continuously updated from interatomic distances, allowing smooth bond breaking and formation. Terms include covalent bonding, valence and torsion angles, lone pairs, hydrogen bonding, Coulomb interactions via charge equilibration, and van der Waals dispersion. ReaxFF is significantly more expensive than simpler potentials (\( \sim 100\times \) LJ cost per atom) but enables simulation of combustion, oxidation, catalysis, and tribochemistry.

Force field selection summary. LJ: noble gases, simple liquids, method benchmarking. EAM: metals and alloys, grain boundaries, dislocations. Tersoff: covalent semiconductors, diamond, amorphous carbon. ReaxFF: reactive chemistry, oxidation, catalysis, wherever bond making/breaking matters. For biological macromolecules: CHARMM, AMBER, OPLS — specialised biomolecular force fields not covered in this course.

Chapter 4: Thermostats, Barostats, and Macroscopic Properties

4.1 Thermodynamic Ensembles in MD

Pure Newtonian MD samples the microcanonical (NVE) ensemble. Real experiments are performed at constant temperature (NVT) or constant temperature and pressure (NPT). Coupling the system to a heat bath or pressure reservoir requires modifying the equations of motion.

NVT (canonical) ensemble. The NVT ensemble fixes particle number \( N \), volume \( V \), and temperature \( T \). The probability of each microstate is proportional to the Boltzmann factor \( e^{-\beta \mathcal{H}} \). An NVT thermostat modifies the equations of motion so that the system samples configurations from this distribution while maintaining dynamical character.

NPT ensemble. The NPT ensemble fixes \( N \), pressure \( P \), and temperature \( T \). The volume fluctuates. An NPT simulation requires both a thermostat and a barostat. The average volume gives the equation of state at the specified \( (T,P) \) conditions; this is essential for computing structural phase transitions and thermal expansion.

4.2 Nosé-Hoover Thermostat

The Nosé-Hoover thermostat introduces an extended Hamiltonian with an additional degree of freedom \( s \) (and its conjugate momentum \( p_s \)) that acts as a friction variable:

\[ \dot{\mathbf{r}}_i = \frac{\mathbf{p}_i}{m_i}, \qquad \dot{\mathbf{p}}_i = \mathbf{F}_i - \xi \mathbf{p}_i, \]\[ \dot{\xi} = \frac{1}{Q}\left(\sum_i \frac{\mathbf{p}_i^2}{m_i} - 3Nk_BT\right), \]

where \( \xi \) is the thermostat friction coefficient and \( Q = 3Nk_BT\tau^2 \) is the thermostat mass (with \( \tau \) the thermostat time constant). When the kinetic energy exceeds \( \frac{3}{2}Nk_BT \), the friction \( \xi \) increases, slowing the particles; when below, it decreases. This feedback loop drives the time-averaged kinetic temperature toward \( T \). The Nosé-Hoover chain (NHC) — introduced by Martyna, Klein, and Tuckerman — chains multiple thermostats to correct ergodicity failures in the single-thermostat case.

The optimal thermostat time constant \( \tau \) is chosen to be comparable to the longest relaxation time in the system — typically 0.1–1 ps for condensed phases.

4.3 Parrinello-Rahman Barostat

The Parrinello-Rahman (PR) barostat couples the simulation cell matrix \( \mathbf{h} = (\mathbf{a}_1 | \mathbf{a}_2 | \mathbf{a}_3) \) to a pressure bath. The cell vectors evolve according to

\[ \ddot{\mathbf{h}} = V\mathbf{W}^{-1} \mathbf{h}^{-T} (\mathbf{\Pi} - P_{\text{ref}}\mathbf{I}), \]

where \( \mathbf{\Pi} \) is the instantaneous stress tensor, \( P_{\text{ref}} \) is the target pressure, and \( \mathbf{W} \) is the barostat mass matrix. The NPT ensemble is realised by combining the PR barostat with the Nosé-Hoover thermostat. This combination allows anisotropic cell deformation, essential for simulating structural phase transitions.

4.4 Radial Distribution Function

The radial distribution function (RDF) \( g(r) \) is one of the most important structural observables extractable from an MD simulation:

\[ g(r) = \frac{V}{N^2} \left\langle \sum_{i \neq j} \delta(r - r_{ij}) \right\rangle. \]

Given a particle at the origin, \( g(r) \) is the ratio of the local number density at distance \( r \) to the bulk number density \( \rho = N/V \). For an ideal gas \( g(r) = 1 \) everywhere. In liquids, oscillatory structure decays to 1 at large \( r \). In crystalline solids, sharp peaks appear at lattice distances.

Example 4.1: Interpreting g(r) for liquid argon.

At \( T = 100 \) K (near the boiling point) and pressure 1 atm, the LJ argon RDF shows:

A hard-core exclusion region with \( g(r) = 0 \) for \( r < \sigma \approx 3.4 \) Å.
A primary coordination shell peak near \( r \approx 3.8 \) Å with \( g(r) \approx 3.0 \), indicating a high probability of nearest neighbours.
A second shell peak near \( r \approx 7 \) Å, reduced and broadened.
Approach to \( g(r) \to 1 \) by \( r \approx 15 \) Å.

The integral \( n_1 = 4\pi\rho \int_0^{r_{\min}} g(r) r^2\, dr \approx 12 \) gives the average coordination number of the first shell, comparable to an FCC crystal. This is consistent with the known tendency of noble-gas liquids near their melting points to retain short-range crystalline order.

4.5 Mean-Square Displacement and Diffusion

Einstein diffusion relation. For a system in equilibrium at sufficiently long times, the mean-square displacement (MSD) of a tagged particle is related to the self-diffusion coefficient \( D \) by

\[ \text{MSD}(t) = \left\langle \left| \mathbf{r}_i(t) - \mathbf{r}_i(0) \right|^2 \right\rangle = 6Dt \quad (t \to \infty), \]

where the factor 6 arises from three independent spatial dimensions. The diffusion coefficient is

\[ D = \lim_{t\to\infty} \frac{\text{MSD}(t)}{6t}. \]

Derivation from the fluctuation-dissipation theorem.

From the Green-Kubo formalism, the diffusion coefficient is related to the velocity autocorrelation function (VACF):

\[ D = \frac{1}{3} \int_0^{\infty} \langle \mathbf{v}_i(t) \cdot \mathbf{v}_i(0) \rangle\, dt. \]

The MSD can be written as

\[ \text{MSD}(t) = \left\langle \left| \int_0^t \mathbf{v}_i(t')\, dt' \right|^2 \right\rangle = 2 \int_0^t (t - t') C_v(t')\, dt', \]

where \( C_v(t') = \langle \mathbf{v}_i(t') \cdot \mathbf{v}_i(0) \rangle \) is the VACF. At long times, if \( C_v(t) \) decays to zero sufficiently fast:

\[ \frac{d(\text{MSD})}{dt} = 2\int_0^t C_v(t')\, dt' \xrightarrow{t \to \infty} 2\int_0^\infty C_v(t')\, dt' = 6D. \]

Integrating: \( \text{MSD}(t) \sim 6Dt \). This connects the macroscopic transport coefficient \( D \) to microscopic velocity fluctuations, exemplifying the fluctuation-dissipation theorem.

Example 4.2: Computing the diffusion coefficient from MSD in LAMMPS.

After equilibration, add to the LAMMPS input:

compute msd all msd com yes
fix msd_out all ave/time 1 1 100 c_msd[4] file msd.dat

The column c_msd[4] is the total MSD (sum of \( x \), \( y \), \( z \) components). Plot MSD versus time and perform a linear regression on the long-time region, avoiding the ballistic regime (early times where \( \text{MSD} \propto t^2 \)). The slope gives \( 6D \). For liquid argon at 100 K, the experimental value is \( D \approx 1.7 \times 10^{-9} \) m\(^2\)/s, readily reproduced by LJ MD to within 10%.

4.6 Pressure and the Virial Theorem

Virial theorem for pressure. For a classical system of \( N \) particles with pairwise interactions, the instantaneous pressure is

\[ P = \frac{Nk_BT}{V} + \frac{1}{3V} \left\langle \sum_{i < j} \mathbf{r}_{ij} \cdot \mathbf{F}_{ij} \right\rangle, \]

where \( \mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j \) and \( \mathbf{F}_{ij} \) is the force on particle \( i \) due to particle \( j \). The first term is the ideal-gas contribution; the second is the virial correction from interactions.

The stress tensor components are

\[ \sigma_{\alpha\beta} = -\frac{1}{V}\left[ \sum_i m_i v_{i\alpha} v_{i\beta} + \sum_{ifrom which the pressure \( P = -\frac{1}{3}\text{tr}(\boldsymbol{\sigma}) \). Elastic constants are computed from stress-strain fluctuations or by applying homogeneous deformations and measuring the stress response.

Example 4.3: Elastic constant calculation via stress-strain in LAMMPS.

To compute the bulk modulus \( B \) of a metal, apply a series of isotropic volume compressions and expansions (LAMMPS deform command). Fit the total energy versus volume to the Birch-Murnaghan equation of state:

\[ E(V) = E_0 + \frac{9V_0 B}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]^3 B' + \left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]^2 \left[6 - 4\left(\frac{V_0}{V}\right)^{2/3}\right] \right\}. \]

Alternatively, compute the stress tensor at each deformation and extract \( C_{ij} \) from Hooke’s law. For FCC copper with the Mishin EAM potential: \( B \approx 138 \) GPa (experiment: 137 GPa), validating the potential.

Chapter 5: Monte Carlo Methods

5.1 Statistical Mechanics Foundations

Monte Carlo simulation is rooted in the observation that thermodynamic ensemble averages,

\[ \langle A \rangle = \frac{\int d\mathbf{r}^N A(\mathbf{r}^N) e^{-\beta U(\mathbf{r}^N)}}{\int d\mathbf{r}^N e^{-\beta U(\mathbf{r}^N)}}, \]

are multi-dimensional integrals that cannot be evaluated analytically. Simple Monte Carlo — drawing configurations uniformly at random from phase space — is catastrophically inefficient: the Boltzmann weight \( e^{-\beta U} \) is vanishingly small for the vast majority of configurations (those with atoms overlapping). Importance sampling, introduced by Metropolis et al. in 1953, solves this by constructing a Markov chain that samples configurations with probability proportional to \( e^{-\beta U(\mathbf{r}^N)} \).

5.2 Markov Chains and Detailed Balance

A Markov chain is a sequence of states \( \mathbf{r}^N_0, \mathbf{r}^N_1, \mathbf{r}^N_2, \ldots \) generated by a stochastic transition rule that depends only on the current state. For the chain to sample the canonical distribution as its stationary distribution, the transition matrix \( \pi(\mathbf{r} \to \mathbf{r}') \) must satisfy detailed balance.

Detailed balance. A Markov chain with stationary distribution \( P(\mathbf{r}) \) satisfies detailed balance if for all pairs of states \( \mathbf{r}, \mathbf{r}' \):

\[ P(\mathbf{r})\, \pi(\mathbf{r} \to \mathbf{r}') = P(\mathbf{r}')\, \pi(\mathbf{r}' \to \mathbf{r}). \]

Detailed balance ensures that the net probability flux between any two states vanishes at equilibrium, guaranteeing that \( P \) is the stationary distribution of the chain.

5.3 The Metropolis Algorithm

Metropolis acceptance criterion. Given the current configuration \( \mathbf{r}^N \) with energy \( U_{\text{old}} \), a trial configuration \( \mathbf{r}'^N \) is generated by randomly displacing one particle. The energy change is \( \Delta U = U_{\text{new}} - U_{\text{old}} \). The trial configuration is accepted with probability

\[ p_{\text{acc}} = \min\!\left(1,\, e^{-\beta \Delta U}\right). \]

If accepted, the new configuration is \( \mathbf{r}'^N \); if rejected, the old configuration is retained and counted again in the ensemble average.

Metropolis satisfies detailed balance.

We must verify \( P(\mathbf{r})\pi(\mathbf{r} \to \mathbf{r}') = P(\mathbf{r}')\pi(\mathbf{r}' \to \mathbf{r}) \) for the canonical distribution \( P(\mathbf{r}) \propto e^{-\beta U(\mathbf{r})} \).

The transition probability decomposes as \( \pi(\mathbf{r} \to \mathbf{r}') = T(\mathbf{r} \to \mathbf{r}') \cdot p_{\text{acc}}(\mathbf{r} \to \mathbf{r}') \), where \( T \) is the probability of proposing the move. For symmetric proposals (\( T(\mathbf{r} \to \mathbf{r}') = T(\mathbf{r}' \to \mathbf{r}) \)), we need only check the acceptance ratio.

Case 1: \( \Delta U \leq 0 \) (move is energetically downhill).

\( p_{\text{acc}}(\mathbf{r} \to \mathbf{r}') = 1 \) and \( p_{\text{acc}}(\mathbf{r}' \to \mathbf{r}) = e^{-\beta(U_{\text{old}} - U_{\text{new}})} = e^{\beta \Delta U} \).

LHS: \( e^{-\beta U_{\text{old}}} \cdot T \cdot 1 \). RHS: \( e^{-\beta U_{\text{new}}} \cdot T \cdot e^{\beta \Delta U} = e^{-\beta U_{\text{new}}} \cdot T \cdot e^{-\beta(U_{\text{old}} - U_{\text{new}})} = e^{-\beta U_{\text{old}}} \cdot T \). Check.

Case 2: \( \Delta U > 0 \) (move is energetically uphill).

\( p_{\text{acc}}(\mathbf{r} \to \mathbf{r}') = e^{-\beta \Delta U} \) and \( p_{\text{acc}}(\mathbf{r}' \to \mathbf{r}) = 1 \).

LHS: \( e^{-\beta U_{\text{old}}} \cdot T \cdot e^{-\beta \Delta U} = e^{-\beta U_{\text{new}}} \cdot T \). RHS: \( e^{-\beta U_{\text{new}}} \cdot T \cdot 1 \). Check.

Both cases satisfy detailed balance. Since the proposal is ergodic, the Metropolis chain converges to the canonical distribution.

5.4 Practical Implementation

A Monte Carlo cycle consists of \( N \) attempted moves, where \( N \) is the number of particles. The displacement magnitude \( \delta_{\max} \) is tuned so that the acceptance ratio is approximately 0.3–0.5. Too large \( \delta_{\max} \) leads to mostly rejections (poor exploration); too small leads to mostly acceptances but tiny moves (also poor exploration). This is tuned during equilibration.

The simulation proceeds in two phases:

Equilibration: \( 10^5 \)–\( 10^6 \) cycles during which the system relaxes to the stationary distribution. Observables accumulated during equilibration are discarded.
Production: Observables are accumulated. Block averaging is used to estimate statistical errors: the production run is divided into \( M \) blocks, the observable is averaged within each block, and the standard deviation of block averages estimates the error.

Example 5.1: MC simulation of a phase transition in a 2D Ising model.

The 2D Ising model on an \( L \times L \) square lattice with spins \( s_i = \pm 1 \) and Hamiltonian \( \mathcal{H} = -J \sum_{\langle ij \rangle} s_i s_j \) undergoes a second-order phase transition at \( T_c = 2J / (k_B \ln(1+\sqrt{2})) \approx 2.269 J/k_B \). A single-spin-flip Metropolis MC simulation with periodic boundary conditions accurately locates this transition. The magnetisation \( m = \langle |N^{-1}\sum_i s_i| \rangle \) serves as the order parameter. Near \( T_c \), finite-size scaling of \( m \) and the susceptibility \( \chi \propto \langle m^2 \rangle - \langle m \rangle^2 \) allows accurate determination of \( T_c \) and critical exponents.

5.5 Applications: Adsorption and Defect Modelling

Monte Carlo naturally handles open-system problems. Adsorption isotherms — the equilibrium uptake of a gas molecule as a function of partial pressure — are among the most important quantities for gas separation membrane and porous material design. Point defect concentrations in crystals at finite temperature are governed by

\[ c_{\text{vac}} = \exp\!\left(-\frac{G_f}{k_B T}\right), \]

where \( G_f = H_f - TS_f \) is the Gibbs free energy of formation. MC allows direct computation of defect equilibria in mixed systems where multiple defect species interact.

5.6 Advanced Sampling

Standard Metropolis MC may fail to efficiently sample configuration space when states are separated by high energy barriers. Two important remedies are:

Umbrella sampling: Introduce a biasing potential \( w(\xi) \) along a reaction coordinate \( \xi \), suppressing high-barrier regions. The unbiased ensemble average is recovered by reweighting:

\[ \langle A \rangle = \frac{\langle A e^{\beta w(\xi)} \rangle_{\text{biased}}}{\langle e^{\beta w(\xi)} \rangle_{\text{biased}}}. \]

Parallel tempering (replica exchange MC): Run \( M \) independent replicas at temperatures \( T_1 < T_2 < \cdots < T_M \). Periodically, attempt to swap configurations between adjacent temperatures with acceptance probability \( \min(1, e^{(\beta_i - \beta_{i+1})(U_i - U_{i+1})}) \). High-temperature replicas cross barriers easily; configurations are passed down to low temperature, allowing the cold replica to escape metastable states.

Chapter 6: Grand Canonical and Kinetic Monte Carlo

6.1 Grand Canonical Monte Carlo

The grand canonical ensemble (\( \mu VT \)) is the natural framework for adsorption, because experiments are performed at fixed chemical potential (partial pressure) rather than fixed particle number.

Chemical potential. The chemical potential \( \mu \) of species \( \alpha \) is the thermodynamic work required to add one particle of species \( \alpha \) to the system at constant \( T \), \( V \), and other particle numbers:

\[ \mu_\alpha = \left(\frac{\partial F}{\partial N_\alpha}\right)_{T,V,N_{\beta \neq \alpha}}. \]

For an ideal gas, \( \mu = \mu^\circ + k_BT\ln(P/P^\circ) \), so fixing \( \mu \) is equivalent to fixing the pressure reservoir.

6.1.1 GCMC Move Types

GCMC Monte Carlo uses three types of moves:

Displacement: Move a randomly selected particle; accept/reject as standard Metropolis.
Insertion: Attempt to insert a particle at a random position. The acceptance probability is

\[ p_{\text{ins}} = \min\!\left(1,\, \frac{V}{N+1} \Lambda^{-3} e^{\beta\mu} e^{-\beta \Delta U}\right), \]

where \( \Lambda = h/\sqrt{2\pi m k_B T} \) is the thermal de Broglie wavelength.

Deletion: Attempt to remove a randomly selected particle. The acceptance probability is

\[ p_{\text{del}} = \min\!\left(1,\, \frac{N}{V} \Lambda^3 e^{-\beta\mu} e^{-\beta \Delta U}\right). \]

Insertion and deletion moves satisfy detailed balance in the grand canonical ensemble. At low pressures (dilute gas), most insertions are rejected due to steric clashes, making GCMC of dense systems challenging. Configurational-bias MC (CBMC) improves efficiency by growing molecules segment by segment, biasing the insertion toward energetically favourable configurations.

Example 6.1: GCMC for CO\(_2\) adsorption in a metal-organic framework (conceptual).

A metal-organic framework (MOF) such as IRMOF-1 has a highly porous structure (\( \sim 80\% \) void fraction) that selectively adsorbs CO\(_2\). A GCMC simulation with the RASPA software package proceeds as follows:

Load the crystal structure of IRMOF-1 as a rigid host framework.
Parametrise CO\(_2\) as a three-site linear molecule with TraPPE force field (\( \varepsilon_C/k_B = 27\) K, \( \sigma_C = 2.8\) Å; quadrupole moment included).
Set up a supercell (\( 2 \times 2 \times 2 \) unit cells) with PBC.
Specify \( T = 298 \) K and a range of fugacities from 0.1 to 50 bar (related to partial pressure by the fugacity coefficient).
Run \( 10^5 \) equilibration cycles followed by \( 10^6 \) production cycles, with equal probability of displacement, insertion, and deletion.
The output is the average number of CO\(_2\) molecules adsorbed, converted to loading in mol/kg, plotted versus pressure to give the adsorption isotherm.

The resulting isotherm is compared against experimental volumetric measurements to validate the force field. This workflow is now standard in virtual screening of MOF libraries for carbon capture.

RASPA software tips. RASPA is open-source and purpose-built for GCMC of porous materials. Key input parameters: ExternalTemperature, ExternalPressure, NumberOfCycles, and the molecule definition files. The framework is typically treated as rigid to reduce cost; flexible-framework GCMC is possible but requires a specialised force field. Always verify that your simulation cell is at least \( 2r_c \) in each periodic direction, with \( r_c \) typically 12–14 Å for Lennard-Jones interactions in MOF simulations.

6.2 Kinetic Monte Carlo

Standard equilibrium MC and MD cannot access timescales longer than \( \sim 1 \) µs for atomic systems. Many technologically important processes — dopant diffusion in semiconductors, grain boundary migration, thin film growth, radiation damage annealing — are dominated by rare thermally activated events with rates far too slow for conventional simulation.

Kinetic Monte Carlo (kMC) bypasses this by treating the dynamics as a continuous-time Markov chain over a discrete set of states, with transition rates given by transition state theory.

6.2.1 Rate Theory and Arrhenius Rates

For a thermally activated process with energy barrier \( E_a \), the rate in the Arrhenius form is

\[ k = \nu_0 \exp\!\left(-\frac{E_a}{k_B T}\right), \]

where \( \nu_0 \) is a prefactor (attempt frequency, typically \( 10^{12} \)–\( 10^{13} \) s\(^{-1}\) for atomic processes). The rates are assumed known — from experiment, DFT nudged elastic band calculations, or transition state theory.

6.2.2 The BKL Algorithm

BKL (Bortz-Kalos-Lebowitz) algorithm correctness. The BKL algorithm generates an exact continuous-time trajectory of a Markov chain with time-independent rates. Given a system with \( M \) possible events and rates \( k_1, k_2, \ldots, k_M \), define the total rate \( R = \sum_{\ell=1}^M k_\ell \). The BKL algorithm proceeds as follows:

Compute \( R = \sum_\ell k_\ell \).
Draw a uniform random number \( u_1 \in (0,1) \). Advance time by \( \Delta t = -\ln(u_1)/R \).
Draw \( u_2 \in (0,1) \). Select event \( j \) such that \( \sum_{\ell=1}^{j-1} k_\ell < u_2 R \leq \sum_{\ell=1}^{j} k_\ell \).
Execute event \( j \). Update the rate list. Return to step 1.

This procedure correctly samples the exponential waiting time and selects events proportional to their rates, producing exact kinetics of the continuous-time master equation.

The time increment \( \Delta t = -\ln(u_1)/R \) follows from the exponential distribution of waiting times: given rate \( R \), the probability that the system remains in the current state for time \( t \) is \( e^{-Rt} \). The logarithm transforms a uniform random number into an exponentially distributed waiting time.

Example 6.2: kMC simulation of surface diffusion.

Consider adatom diffusion on an FCC (100) surface. An adatom at site \( i \) can hop to any of its four nearest-neighbour empty sites. The rate for each hop depends on the local coordination:

\[ k_{\text{hop}} = \nu_0 \exp\!\left(-\frac{E_a + n_n \cdot \Delta E}{k_B T}\right), \]

where \( n_n \) is the number of lateral neighbours and \( \Delta E \approx 0.1\text{–}0.2 \) eV penalises breaking lateral bonds. At \( T = 500 \) K with \( \nu_0 = 10^{13} \) s\(^{-1}\), \( E_a = 0.5 \) eV, the typical hop rate is \( \sim 10^8 \) s\(^{-1}\), corresponding to a mean waiting time of 10 ns — inaccessible to classical MD. The kMC simulation reaches seconds of physical time easily, allowing computation of tracer diffusion coefficients and island nucleation statistics in thin film growth.

6.2.3 SPPARKS for Mesoscale kMC

SPPARKS (Stochastic Parallel PARticle Kinetic Simulator, Sandia National Labs) implements on-lattice kMC for grain growth, sintering, irradiation, and deposition problems on parallel architectures. An application model is defined by specifying lattice geometry, site variables, and event rates. SPPARKS uses the BKL algorithm with efficient partial-sum trees for rate selection, scaling to billions of sites.

kMC limitations. kMC requires a complete and closed catalogue of all possible events and their rates. For complex systems (amorphous materials, multi-component alloys), constructing this catalogue is non-trivial. Adaptive kMC methods (ART nouveau, dimer method) autonomously identify transition states, but at significantly greater computational cost. When rates depend on unknown configurations, off-lattice kMC coupled with ML potentials is an active research frontier.

Chapter 7: Density Functional Theory — Theory

7.1 The Many-Body Electronic Structure Problem

The properties of materials at the atomic scale are determined by the behaviour of electrons. The exact many-body Schrödinger equation for \( N_e \) electrons and \( N_I \) nuclei is

\[ \hat{H}\Psi = E\Psi, \qquad \hat{H} = -\sum_i \frac{\hbar^2}{2m_e}\nabla_i^2 + \sum_I \frac{P_I^2}{2M_I} + V_{ee} + V_{eI} + V_{II}, \]

where \( V_{ee} = \frac{e^2}{4\pi\varepsilon_0}\sum_{i

7.1.1 Born-Oppenheimer Approximation

The nuclei are \( 10^3 \)–\( 10^5 \) times heavier than electrons and move orders of magnitude more slowly. To an excellent approximation, the electrons instantaneously adjust to any nuclear configuration. This allows separation of the wavefunction:

\[ \Psi(\mathbf{r}^{N_e}, \mathbf{R}^{N_I}) \approx \psi_e(\mathbf{r}^{N_e}; \mathbf{R}^{N_I}) \cdot \chi(\mathbf{R}^{N_I}). \]

The electronic Schrödinger equation with clamped nuclei is

\[ \hat{H}_e \psi_e = E_e(\mathbf{R}^{N_I}) \psi_e, \qquad \hat{H}_e = -\sum_i \frac{\hbar^2}{2m_e}\nabla_i^2 + V_{ee} + V_{eI}(\mathbf{R}^{N_I}). \]

The electronic energy \( E_e(\mathbf{R}) \) constitutes the potential energy surface for nuclear motion — this is what MD force fields attempt to approximate.

7.1.2 Hartree Approximation

The Hartree approximation replaces the correlated many-electron wavefunction with a product of single-particle orbitals: \( \psi_e \approx \prod_i \phi_i(\mathbf{r}_i) \). This reduces the \( 3N_e \)-dimensional problem to \( N_e \) coupled 3D equations. The Hartree-Fock method further antisymmetrises this product (Slater determinant) to satisfy Pauli exclusion, adding an exact but non-local exchange term. DFT provides a formally exact alternative.

7.2 Hohenberg-Kohn Theorems

The foundation of DFT rests on two theorems proved by Hohenberg and Kohn in 1964.

Hohenberg-Kohn existence theorem (HK1). The external potential \( v_{\text{ext}}(\mathbf{r}) \) is determined uniquely (up to a constant) by the ground-state electron density \( n_0(\mathbf{r}) \). Consequently, since \( v_{\text{ext}} \) determines the Hamiltonian which determines all properties of the system, the ground-state electron density \( n_0(\mathbf{r}) \) determines all ground-state properties.

Proof by contradiction (reductio ad absurdum).

Suppose two different external potentials \( v \) and \( v' = v + \Delta v \) (with \( \Delta v \) not a constant) yield the same ground-state density \( n_0(\mathbf{r}) \), with Hamiltonians \( \hat{H} \) and \( \hat{H}' \) and ground states \( |\psi\rangle \) and \( |\psi'\rangle \).

By the variational principle, \( \langle \psi | \hat{H} | \psi \rangle < \langle \psi' | \hat{H} | \psi' \rangle \), which gives:

\[ E_0 < \langle \psi' | \hat{H} | \psi' \rangle = \langle \psi' | \hat{H}' | \psi' \rangle + \langle \psi' | \hat{H} - \hat{H}' | \psi' \rangle = E_0' + \int n_0(\mathbf{r})\left[v(\mathbf{r}) - v'(\mathbf{r})\right] d\mathbf{r}. \]

By the same argument with \( \psi \) and \( \psi' \) swapped:

\[ E_0' < E_0 - \int n_0(\mathbf{r})\left[v(\mathbf{r}) - v'(\mathbf{r})\right] d\mathbf{r}. \]

Adding these two inequalities yields \( E_0 + E_0' < E_0 + E_0' \), a contradiction. Hence two different potentials cannot produce the same ground-state density, proving HK1.

Hohenberg-Kohn variational theorem (HK2). The ground-state energy \( E_0 \) can be obtained variationally: for any trial density \( \tilde{n}(\mathbf{r}) \geq 0 \) with \( \int \tilde{n}\, d\mathbf{r} = N_e \),

\[ E_0 \leq E[\tilde{n}] = F_{\text{HK}}[\tilde{n}] + \int \tilde{n}(\mathbf{r}) v_{\text{ext}}(\mathbf{r})\, d\mathbf{r}, \]

where \( F_{\text{HK}}[n] = \langle \psi[n] | \hat{T} + \hat{V}_{ee} | \psi[n] \rangle \) is the universal Hohenberg-Kohn functional, independent of the external potential.

The HK theorems establish that the density is the fundamental variable, but provide no practical prescription for computing \( F_{\text{HK}} \). Kohn and Sham resolved this in 1965.

7.3 Kohn-Sham Equations

Kohn-Sham mapping theorem. There exists a unique auxiliary non-interacting system — the Kohn-Sham system — whose ground-state density equals the exact ground-state density of the interacting system. The Kohn-Sham orbitals \( \{\phi_i(\mathbf{r})\} \) and eigenvalues \( \{\varepsilon_i\} \) satisfy the Kohn-Sham equations:

\[ \left[ -\frac{\hbar^2}{2m_e}\nabla^2 + v_{\text{eff}}(\mathbf{r}) \right] \phi_i(\mathbf{r}) = \varepsilon_i \phi_i(\mathbf{r}), \]

where the effective potential is

\[ v_{\text{eff}}(\mathbf{r}) = v_{\text{ext}}(\mathbf{r}) + v_H(\mathbf{r}) + v_{xc}(\mathbf{r}), \]

with the Hartree potential \( v_H(\mathbf{r}) = e^2 \int \frac{n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}\, d\mathbf{r}' \) and the exchange-correlation potential \( v_{xc}(\mathbf{r}) = \delta E_{xc}[n] / \delta n(\mathbf{r}) \). The ground-state density is \( n(\mathbf{r}) = \sum_i^{N_e} |\phi_i(\mathbf{r})|^2 \).

Kohn-Sham equations. The Kohn-Sham equations are a set of \( N_e \) coupled single-particle Schrödinger equations that must be solved self-consistently, because \( v_{\text{eff}} \) depends on \( n(\mathbf{r}) \), which in turn depends on the orbitals \( \phi_i \). The exchange-correlation energy \( E_{xc}[n] \) captures all non-classical electron-electron interactions not included in the Hartree term.

Kohn-Sham equations from variational principle.

The total energy functional is

\[ E[n] = T_s[n] + E_H[n] + E_{xc}[n] + \int n(\mathbf{r}) v_{\text{ext}}(\mathbf{r})\, d\mathbf{r}, \]

where \( T_s[n] = -\frac{\hbar^2}{2m_e}\sum_i \int \phi_i^* \nabla^2 \phi_i\, d\mathbf{r} \) is the kinetic energy of the non-interacting Kohn-Sham system and \( E_H[n] = \frac{e^2}{2}\int\int \frac{n(\mathbf{r})n(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\, d\mathbf{r}\, d\mathbf{r}' \) is the Hartree energy.

Minimising \( E[n] \) with respect to \( \phi_i^* \) subject to orthonormality constraints \( \langle \phi_i | \phi_j \rangle = \delta_{ij} \) via Lagrange multipliers gives:

\[ \frac{\delta}{\delta \phi_i^*}\left[ E - \sum_{j} \varepsilon_j \langle \phi_j | \phi_j \rangle \right] = 0 \]\[ \implies -\frac{\hbar^2}{2m_e}\nabla^2 \phi_i + \left[ v_{\text{ext}}(\mathbf{r}) + v_H(\mathbf{r}) + v_{xc}(\mathbf{r}) \right]\phi_i = \varepsilon_i \phi_i, \]

identifying \( v_H = \delta E_H / \delta n \) and \( v_{xc} = \delta E_{xc} / \delta n \). This is the Kohn-Sham equation. The formalism is exact if the true \( E_{xc}[n] \) were known; the approximation enters only through the exchange-correlation functional.

7.4 Self-Consistent Field Procedure

The Kohn-Sham equations must be solved iteratively:

Guess an initial electron density \( n^{(0)}(\mathbf{r}) \) (e.g., superposition of atomic densities).
Compute \( v_H \) and \( v_{xc} \), construct \( v_{\text{eff}} \).
Solve the KS eigenvalue problem to obtain new orbitals \( \{\phi_i\}^{(k)} \).
Compute new density \( n^{(k)}(\mathbf{r}) = \sum_i |\phi_i^{(k)}|^2 \).
Mix old and new density (Pulay or Anderson mixing in practice).
Check convergence: \( \|n^{(k+1)} - n^{(k)}\| < \varepsilon_{\text{conv}} \). If not converged, return to step 2.

Convergence tolerances: total energy \( < 10^{-6} \) Ry per atom; forces \( < 10^{-3} \) Ry/bohr for structural relaxation.

7.5 Exchange-Correlation Functionals

Exchange-correlation functional. The exchange-correlation functional \( E_{xc}[n] \) accounts for (1) quantum mechanical exchange energy arising from the antisymmetry of the wavefunction, and (2) correlation energy from the correlated motion of electrons. Exact expressions exist only for the homogeneous electron gas; all practical functionals involve controlled approximations.

Local Density Approximation (LDA): The exchange-correlation energy density at each point is taken from the homogeneous electron gas at the same density:

\[ E_{xc}^{\text{LDA}}[n] = \int n(\mathbf{r}) \varepsilon_{xc}^{\text{HEG}}(n(\mathbf{r}))\, d\mathbf{r}. \]

LDA systematically overbinds (overestimates cohesive energies by \( \sim 10\% \)) and underestimates band gaps by 30–50%.

Generalised Gradient Approximation (GGA-PBE): Adds dependence on the gradient of the density \( |\nabla n| \):

\[ E_{xc}^{\text{GGA}}[n] = \int n(\mathbf{r}) \varepsilon_{xc}^{\text{GGA}}\!\left(n(\mathbf{r}), |\nabla n(\mathbf{r})|\right) d\mathbf{r}. \]

The PBE functional (Perdew, Burke, Ernzerhof, 1996) is the most widely used GGA. Band gaps remain underestimated.

Hybrid Functionals (HSE06): Introduce a fraction \( \alpha \) of exact (Hartree-Fock) exchange:

\[ E_{xc}^{\text{hyb}} = \alpha E_x^{\text{HF}} + (1-\alpha) E_x^{\text{GGA}} + E_c^{\text{GGA}}. \]

HSE06 (Heyd, Scuseria, Ernzerhof) uses a screened Coulomb kernel making it tractable for periodic systems. Band gaps are typically accurate to \( \pm 0.2 \) eV for semiconductors.

Functional choice guide. For geometry optimisation, phonons, and surface energies: PBE. For electronic structure and band gaps: HSE06 or GW. For strongly correlated systems (transition metal oxides, f-electron compounds): DFT+U with Hubbard \( U \) correction. Never use LDA/GGA band gaps as quantitative predictions without noting the systematic underestimation.

7.6 Basis Sets and Pseudopotentials

7.6.1 Plane-Wave Basis

Bloch’s theorem states that wavefunctions of a periodic crystal are Bloch functions: \( \phi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}} u_{n\mathbf{k}}(\mathbf{r}) \), where \( u_{n\mathbf{k}} \) is periodic. Expanding in a plane-wave basis:

\[ u_{n\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{G}} c_{n\mathbf{k}}^{\mathbf{G}} e^{i\mathbf{G}\cdot\mathbf{r}}, \]

truncated at kinetic energy cutoff \( E_{\text{cut}} = \frac{\hbar^2}{2m_e}|\mathbf{k}+\mathbf{G}|^2 \). Plane waves are systematic (convergence controlled by one parameter), unbiased, and straightforwardly implement PBC. Their disadvantage is the inability to efficiently represent rapidly oscillating wavefunctions near atomic cores.

7.6.2 Pseudopotentials and the PAW Method

Core electrons do not participate significantly in bonding but their rapid oscillations near the nucleus require a large plane-wave cutoff. Pseudopotentials replace the true electron-ion interaction with a smoother effective potential that reproduces scattering properties outside a core radius \( r_c \).

The Projector Augmented Wave (PAW) method (Blöchl, 1994) defines a linear transformation between smooth pseudo-wavefunctions and the true all-electron wavefunctions, allowing recovery of all-electron quantities at a cost similar to ultrasoft pseudopotentials. PAW is the default in VASP and Quantum ESPRESSO.

7.7 k-point Sampling

In a crystalline solid, the Brillouin zone (BZ) integrals must be discretised. The Monkhorst-Pack (MP) scheme generates a uniform grid of \( N_1 \times N_2 \times N_3 \) k-points in the BZ. For metals, finer grids and Methfessel-Paxton or Marzari-Vanderbilt smearing of the Fermi surface are required.

k-point convergence testing. Always verify that your total energy is converged with respect to k-point density. For a semiconductor like silicon, a \( 4\times4\times4 \) MP grid typically converges the total energy to \( < 1 \) meV/atom. For metals, \( 16\times16\times16 \) or denser grids may be needed. Band structure calculations use a fine path through high-symmetry k-points, but this does not need to be self-consistent — use the charge density from the SCF calculation.

Chapter 8: DFT in Practice — Band Structure, DOS, and Properties

8.1 The Quantum ESPRESSO Workflow

Quantum ESPRESSO (QE) is an open-source package for DFT calculations using plane-wave basis sets and pseudopotentials. The primary programs are:

pw.x: self-consistent field calculation and structural relaxation.
dos.x: computes the density of states from a non-self-consistent calculation on a fine k-point grid.
bands.x: unfolds and plots the band structure along specified k-paths.
pp.x: post-processing (charge density, wave functions, electrostatic potential).
ph.x: phonon calculations via density functional perturbation theory.

Example 8.1: DFT convergence testing for silicon.

Silicon (diamond cubic, \( a = 5.431 \) Å) is the canonical DFT benchmark. The SCF input file specifies the cell parameter celldm(1) = 10.263 (in bohr), two atoms per unit cell in the diamond structure, a PBE pseudopotential, and an \( 8\times8\times8 \) Monkhorst-Pack k-grid.

Convergence testing protocol:

Fix k-points at \( 8\times8\times8 \); vary ecutwfc (wavefunction cutoff) from 20 to 80 Ry. Plot total energy versus cutoff; find the value where \( \Delta E < 1 \) meV/atom (typically 40–50 Ry for norm-conserving pseudopotentials, 25–30 Ry for PAW).
Fix ecutwfc at the converged value; vary k-point mesh from \( 2\times2\times2 \) to \( 12\times12\times12 \). Find the mesh where \( \Delta E < 1 \) meV/atom (typically \( 6\times6\times6 \) for Si).
Verify the lattice constant by relaxing the unit cell at the converged parameters. The PBE value \( a \approx 5.47 \) Å overestimates the experimental value (5.43 Å) by \( \sim 0.7\% \), consistent with typical GGA overbinding.

8.2 Band Structure Calculations

Example 8.2: Band structure of silicon (indirect gap).

Silicon has an indirect band gap: the valence band maximum (VBM) is at the \( \Gamma \) point and the conduction band minimum (CBM) is at approximately \( 0.85 \times \overrightarrow{\Gamma X} \) along the \( \Delta \) direction. The workflow:

SCF calculation on an \( 8\times8\times8 \) k-grid to obtain the ground-state charge density.
NSCF (non-self-consistent) calculation using the converged charge density, with k-points along the high-symmetry path \( L \to \Gamma \to X \to U,K \to \Gamma \).
Run bands.x to extract and sort the eigenvalues; post-process with plotband.x.

The PBE band gap of Si is \( \sim 0.6 \) eV (experimental: 1.17 eV at 0 K), illustrating the systematic DFT band gap underestimation. HSE06 gives \( \sim 1.1 \) eV. The band structure reveals the six equivalent conduction band minima related by cubic symmetry, giving silicon its sixfold valley degeneracy exploited in strain engineering and valleytronics.

8.3 Density of States

The electronic density of states is

\[ g(E) = \sum_{n} \int_{\text{BZ}} \frac{d\mathbf{k}}{(2\pi)^3} \delta(E - \varepsilon_{n\mathbf{k}}), \]

computed numerically by replacing the delta function with a Gaussian or Lorentzian broadening \( \sigma = 0.02\text{–}0.1 \) eV. The projected (partial) DOS (PDOS) decomposes \( g(E) \) onto atomic orbitals and angular momentum channels:

\[ g_{\ell, I}(E) = \sum_{n,\mathbf{k}} w_\mathbf{k} |\langle \phi_{\ell,I} | \psi_{n\mathbf{k}} \rangle|^2 \delta(E - \varepsilon_{n\mathbf{k}}), \]

where \( \phi_{\ell,I} \) is an atomic orbital of atom \( I \) with angular momentum \( \ell \). PDOS reveals the character of bands: in TiO\(_2\), the valence band is dominated by O-2p states while the conduction band is Ti-3d.

8.4 Band Gap Engineering

The band gap of a semiconductor can be tuned by:

Biaxial strain: Strain modifies the crystal field and lattice constant. In Si, 1% biaxial tensile strain reduces the band gap by \( \sim 0.1 \) eV and lifts the conduction band valley degeneracy, enhancing electron mobility.
Alloying: In \( \text{In}_{1-x}\text{Ga}_x\text{N} \) alloys, the band gap varies from 0.7 eV (InN) to 3.4 eV (GaN), covering the entire visible spectrum. DFT calculations with special quasi-random structures (SQS) model the alloy disorder.
Quantum confinement: In quantum wells and quantum dots, confinement of electronic wavefunctions blueshifts the effective band gap. DFT slab calculations model the electronic structure of ultrathin films.

8.5 Optical Properties

The imaginary part of the dielectric function \( \varepsilon_2(\omega) \) governs optical absorption:

\[ \varepsilon_2(\omega) = \frac{4\pi^2 e^2}{m_e^2 \omega^2} \sum_{n,m,\mathbf{k}} w_\mathbf{k} |\langle \psi_{n\mathbf{k}} | \hat{\mathbf{e}} \cdot \mathbf{p} | \psi_{m\mathbf{k}} \rangle|^2 \delta(\varepsilon_{m\mathbf{k}} - \varepsilon_{n\mathbf{k}} - \hbar\omega), \]

summing over all vertical transitions from occupied state \( n \) to unoccupied state \( m \) with photon energy \( \hbar\omega \). The real part \( \varepsilon_1(\omega) \) follows from the Kramers-Kronig relations. From \( \varepsilon_{1,2} \), one computes the absorption coefficient, reflectivity, and refractive index. The joint density of states (JDOS) — the convolution of occupied and unoccupied DOS — determines the overall shape of the absorption spectrum.

8.6 Charge Density Analysis

8.6.1 Bader Analysis

Example 8.3: Bader charge analysis of TiO\(_2\).

Bader analysis partitions the electron density among atomic basins using zero-flux surfaces (\( \nabla n \cdot \hat{n} = 0 \) on the surface). The Bader charge on atom \( I \) is

\[ q_I = Z_I - \int_{\Omega_I} n(\mathbf{r})\, d\mathbf{r}, \]

where \( Z_I \) is the nuclear charge and the integral is over the Bader basin \( \Omega_I \). For rutile TiO\(_2\) with PBE: Ti acquires \( +2.4e \) and each O acquires \( -1.2e \) (formal oxidation states Ti\(^{4+}\), O\(^{2-}\), but DFT gives fractional values). Bader analysis identifies charge transfer at interfaces and surface sites with unusual oxidation states, guiding catalyst design. The Henkelman group code computes Bader volumes from the QE charge density output in CUBE format.

8.6.2 Electron Localisation Function

The electron localisation function (ELF) ranges from 0 to 1 and indicates regions of localised electron pairs. ELF \( \approx 1 \) at covalent bonds and lone pairs; ELF \( \approx 0.5 \) in electron-gas-like regions; ELF \( \approx 0 \) at shell boundaries. ELF is visualised on a colour scale with VESTA or XCrysDen.

Quantum ESPRESSO workflow tips. Always use the same pseudopotential library consistently within a study. The SSSP (Standard Solid State Pseudopotential) library is recommended for materials benchmarking. For charged systems (defects, surfaces), use a compensating background charge and apply the Freysoldt finite-size correction to the defect formation energy. For metallic systems, the smearing keyword and choice of smearing method (Methfessel-Paxton recommended) strongly affect convergence — do not use Gaussian smearing width larger than 0.05 Ry for production calculations.

Chapter 9: Machine Learning for Materials Science

9.1 Motivation: The Accuracy-Cost Dilemma

Classical interatomic potentials — LJ, EAM, Tersoff — are computationally cheap (\( \sim 10^{-4} \) s/atom/step) and enable microsecond simulations of million-atom systems. However, they are non-transferable: a potential parameterised for FCC copper may fail catastrophically at grain boundaries, surfaces, or in alloy configurations far from the training set. DFT provides accurate energies and forces for essentially any configuration but at enormous cost (\( O(N_e^3) \), \( \sim 10^4 \) atom·step/second for 100 atoms). The gap spans eight orders of magnitude.

Machine-learning interatomic potentials (MLIPs) learn the DFT potential energy surface from a training set of DFT-labelled configurations. Inference cost is \( \sim 10^{-2} \)–\( 10^{0} \) s/atom/step — two to four orders of magnitude faster than DFT — while maintaining DFT-level accuracy for configurations within the training distribution.

9.2 Atom-Centred Symmetry Functions: Behler-Parrinello NNPs

The earliest systematic MLIP framework, introduced by Behler and Parrinello (2007), represents the total energy as a sum of atomic contributions:

\[ E = \sum_i E_i(\mathbf{G}_i), \]

where \( \mathbf{G}_i \) is a vector of atom-centred symmetry functions (ACSFs) encoding the local chemical environment of atom \( i \). ACSFs are constructed to be invariant under translation, rotation, and permutation of identical atoms.

Radial functions:

\[ G_i^2 = \sum_{j \neq i} e^{-\eta (r_{ij} - r_s)^2} f_c(r_{ij}), \]

Angular functions:

\[ G_i^4 = 2^{1-\zeta} \sum_{j \neq i} \sum_{k \neq i,j} (1 + \lambda \cos\theta_{ijk})^\zeta e^{-\eta(r_{ij}^2 + r_{ik}^2 + r_{jk}^2)} f_c(r_{ij}) f_c(r_{ik}) f_c(r_{jk}), \]

with cutoff function \( f_c(r) = \frac{1}{2}[\cos(\pi r/r_c) + 1] \) for \( r \leq r_c \). The vector \( \mathbf{G}_i \) is fed into an element-specific multilayer perceptron (MLP) to predict \( E_i \). Forces are obtained analytically as \( \mathbf{F}_i = -\partial E / \partial \mathbf{r}_i \) via chain rule through the ACSF definition.

9.3 Graph Neural Networks for Materials

The Behler-Parrinello approach requires hand-crafted symmetry functions. Graph neural networks (GNNs) learn feature representations automatically from the graph structure of the material.

Graph neural network. A graph neural network is a class of neural network that operates on graph-structured data. For a material represented as a graph \( G = (V, E) \) where nodes \( V \) are atoms and edges \( E \) connect atom pairs within a cutoff radius, a GNN learns a mapping from node and edge features to a target property (energy, band gap, etc.) by iteratively aggregating information from the local neighbourhood of each node.

9.3.1 Graph Representation of Crystal Structures

A crystal or molecule is represented as a graph:

Nodes: each atom \( i \) has an initial feature vector \( \mathbf{h}_i^{(0)} \) encoding atomic number (one-hot or embedding), oxidation state, electronegativity.
Edges: atom pairs \( (i,j) \) with \( r_{ij} < r_c \) are connected. Edge features \( \mathbf{e}_{ij} \) encode distance, direction (unit vector), and potentially bond type.
For periodic crystals, edges include periodic images, and edge vectors carry the full geometric information \( \mathbf{r}_{ij} = \mathbf{r}_j + \mathbf{L}\mathbf{n} - \mathbf{r}_i \) where \( \mathbf{n} \) is the image vector.

9.3.2 Message Passing Framework

Message passing. The message-passing neural network (MPNN) framework updates node features over \( L \) rounds (layers). At each layer \( \ell \):

Message computation:

\[ \mathbf{m}_{ij}^{(\ell)} = \phi_{\text{msg}}\!\left(\mathbf{h}_i^{(\ell)}, \mathbf{h}_j^{(\ell)}, \mathbf{e}_{ij}\right) \]

Aggregation:

\[ \mathbf{m}_i^{(\ell)} = \sum_{j \in \mathcal{N}(i)} \mathbf{m}_{ij}^{(\ell)} \]

Node update:

\[ \mathbf{h}_i^{(\ell+1)} = \phi_{\text{upd}}\!\left(\mathbf{h}_i^{(\ell)}, \mathbf{m}_i^{(\ell)}\right) \]

After \( L \) layers, a readout function \( \mathbf{R} \) pools node features to predict the target:

\[ \hat{y} = \mathbf{R}\!\left(\left\{ \mathbf{h}_i^{(L)} \right\}_{i \in V}\right) = \sum_i \phi_{\text{out}}\!\left(\mathbf{h}_i^{(L)}\right). \]

Equivariance. A function \( f: X \to Y \) between spaces with group \( G \) actions is G-equivariant if \( f(g \cdot x) = g \cdot f(x) \) for all \( g \in G \). For a physical force field, the potential energy must be invariant (\( E \) is a scalar: \( f(gx) = f(x) \)) while the forces must be equivariant under rotation (\( \mathbf{F}(R\mathbf{r}^N) = R\mathbf{F}(\mathbf{r}^N) \)). Architectures that build in E(3) equivariance — the symmetry group of 3D Euclidean space including rotations, reflections, and translations — are physically correct by construction and achieve dramatically improved data efficiency.

9.3.3 Invariant GNNs: SchNet

SchNet (Schütt et al., 2018) achieves rotation invariance by using only scalar interatomic distances as edge features. Continuous-filter convolutional layers compute:

\[ \mathbf{h}_i^{(\ell+1)} = \sum_{j \in \mathcal{N}(i)} \mathbf{h}_j^{(\ell)} \odot W\!\left(r_{ij}\right), \]

where \( W(r_{ij}) = \phi_\theta(\mathbf{e}(r_{ij})) \) is a distance-dependent filter network applied element-wise. The radial basis expansion \( \mathbf{e}(r) \) uses Gaussian functions centred at a range of distances. SchNet is invariant to rotation by construction but cannot distinguish mirror-image configurations, limiting its accuracy for chiral systems.

9.4 Equivariant Networks: NequIP and MACE

Equivariant GNNs explicitly transform geometric quantities (vectors, tensors) under rotation, enabling learning of directional interactions.

NequIP (Neural Equivariant Interatomic Potentials, Batzner et al., 2022) uses irreducible representations of SO(3) — spherical harmonics — as features. Node features at each layer are tensors of defined rotational degree \( \ell \): scalars (\( \ell = 0 \)), vectors (\( \ell = 1 \)), and higher-order tensors. Tensor products (Clebsch-Gordan products) mix features of different orders while maintaining equivariance. NequIP achieves DFT accuracy with \( 10\text{–}100\times \) fewer training configurations than non-equivariant networks.

MACE (Multi-ACE, Batatia et al., 2022) generalises the Atomic Cluster Expansion (ACE) body-order formalism to a message-passing framework. MACE features are equivariant and systematically expandable in body order, achieving state-of-the-art accuracy on benchmark datasets (MD17, Materials Project) at competitive inference speed.

Universal approximation for GNNs on graphs. A GNN with sufficient depth and width, using a sufficiently expressive aggregation function, is a universal approximator of functions on graphs that are invariant to node permutation. For materials, this implies that a suitably designed GNN can approximate any smooth, permutation-invariant potential energy surface to arbitrary accuracy, given sufficient training data and sufficient model capacity.

9.5 Automatic Differentiation and Forces

A key advantage of ML potentials is that forces are obtained analytically from the energy prediction via automatic differentiation (AD), requiring no explicit force formula:

\[ \mathbf{F}_i = -\nabla_{\mathbf{r}_i} E_\theta(\mathbf{r}^N), \]

where \( E_\theta \) is the GNN energy model parametrised by \( \theta \). Reverse-mode AD (backpropagation) computes the gradient of a scalar with respect to all inputs in a single backward pass at cost proportional to the forward pass.

Example 9.1: Automatic differentiation for forces in PyTorch.

The following pseudocode sketch illustrates the essential pattern used in ML potential implementations:

import torch

# positions: (N_atoms, 3), requires_grad=True for AD
pos = torch.tensor(coordinates, dtype=torch.float64,
                   requires_grad=True)

# Forward pass: compute energy from GNN
energy = model(atomic_numbers, pos, edge_index, edge_vectors)

# Backward pass: compute forces via AD
energy_sum = energy.sum()
energy_sum.backward()
forces = -pos.grad  # F_i = -dE/dr_i

For training with a force-matching loss, create_graph=True must be passed to retain the computation graph through the gradient operation, enabling differentiation of forces with respect to model parameters.

9.6 Training Workflow

9.6.1 Dataset Curation

A high-quality training dataset must:

Cover the configuration space of the intended application: include equilibrium structures, strained cells, defects, surfaces, high-temperature snapshots from ab initio MD (AIMD), and transition-state configurations.
Be DFT-consistent: use the same functional (typically PBE or r\(^2\)SCAN), pseudopotential, k-point grid, and convergence settings throughout.
Be size-representative: train on 50–500-atom cells to capture medium-range interactions.

Active learning efficiently extends the dataset: run MLIP MD, identify configurations where model uncertainty is high, add them to the training set, retrain. Iteration continues until uncertainty is acceptable.

9.6.2 Loss Function

The training loss combines energy and force targets:

\[ \mathcal{L}(\theta) = \lambda_E \sum_s \frac{(E_\theta^{(s)} - E_{\text{DFT}}^{(s)})^2}{N_s^2} + \lambda_F \sum_s \frac{1}{3N_s}\sum_i \|\mathbf{F}_{\theta,i}^{(s)} - \mathbf{F}_{\text{DFT},i}^{(s)}\|^2, \]

where the sum is over training structures \( s \), \( N_s \) is the number of atoms in structure \( s \), and \( \lambda_E, \lambda_F \) are weighting hyperparameters. Forces contain \( 3N \) data points per structure versus one for energy; the weighting typically favours forces (\( \lambda_F \gg \lambda_E \)). Stress tensor targets can be included for NPT MD.

9.6.3 Validation and Uncertainty Quantification

Training, validation, and test splits should be stratified: do not put configurations from the same AIMD trajectory into both training and test sets. Target validation metrics:

Energy RMSE: \( < 1 \) meV/atom.
Force RMSE: \( < 50 \) meV/Å.

Uncertainty quantification (UQ) estimates model confidence for out-of-distribution configurations. Common approaches: deep ensembles (train \( M \) independent models; standard deviation of predictions is the uncertainty), Monte Carlo dropout, or last-layer variance in Gaussian process regression.

Example 9.2: ML potential training loop pseudocode.

The following pseudocode sketches the standard training loop for an MLIP in PyTorch. The model is trained to minimize a combined energy-force loss; a learning rate scheduler reduces the step size when validation loss plateaus, a common strategy for training deep materials models.

model = MACE(max_ell=2, num_channels=128, num_layers=4)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer)

for epoch in range(max_epochs):
    model.train()
    for batch in train_loader:
        optimizer.zero_grad()

        # Forward pass
        E_pred = model(batch)
        # Forces via reverse-mode AD
        F_pred = compute_forces(E_pred, batch.pos)

        # Loss function
        dE = (E_pred / batch.natoms - batch.energy / batch.natoms)
        dF = (F_pred - batch.forces)
        loss_E = (dE ** 2).mean()
        loss_F = (dF ** 2).mean()
        loss = lambda_E * loss_E + lambda_F * loss_F

        # Backward pass and gradient step
        loss.backward()
        torch.nn.utils.clip_grad_norm_(model.parameters(), 10.0)
        optimizer.step()

    # Evaluate on validation set
    model.train(False)
    val_loss = compute_validation_loss(model, val_loader)
    scheduler.step(val_loss)

    if val_loss < best_val_loss:
        best_val_loss = val_loss
        torch.save(model.state_dict(), 'best_model.pt')

9.7 GNN Forward Pass: Worked Example

Example 9.3: GNN forward pass on a small molecule graph (water dimer).

Consider the water dimer: two H\(_2\)O molecules, 6 atoms (O\(_1\), H\(_2\), H\(_3\), O\(_4\), H\(_5\), H\(_6\)). Build the graph with cutoff \( r_c = 4 \) Å.

Step 1: Initialise node features. Each atom gets an embedding: O maps to a 64-dimensional learned vector, H maps to a different 64-dimensional vector, looked up from an embedding table indexed by atomic number.

Step 2: Build edge list. All pairs within 4 Å are connected: intra-molecular O-H bonds (both molecules), inter-molecular O-O, O-H, H-H contacts if within cutoff. Compute edge vectors \( \mathbf{e}_{ij} = \mathbf{r}_j - \mathbf{r}_i \) and distances \( r_{ij} \).

Step 3: Edge features. Expand \( r_{ij} \) in a radial basis:

\[ e_{ij,k} = \exp\!\left(-\eta_k(r_{ij} - \mu_k)^2\right), \quad k = 1,\ldots,20, \]

a 20-dimensional Gaussian basis with centres spanning \( [0, r_c] \).

Step 4: Message passing (layer 1). For atom O\(_1\): collect messages from neighbours H\(_2\), H\(_3\), and any inter-molecular contacts within cutoff. Aggregate by summing. Update O\(_1\)’s feature vector using a gated network or residual MLP.

Step 5: Readout. After \( L = 3 \) layers, apply a two-layer MLP to each atom’s feature vector to produce atomic energy contribution \( E_i \). Total energy: \( E = \sum_i E_i \).

For an equivariant model (NequIP/MACE), Step 3 uses spherical harmonic projections of \( \hat{\mathbf{r}}_{ij} \) as directional features, and the MLP is replaced by Clebsch-Gordan tensor product operations that preserve SO(3) equivariance, enabling the model to learn the orientation of the hydrogen bond between the two water molecules.

9.8 Property Prediction Beyond Forces

The same GNN framework applies to predicting any material property that depends on atomic structure:

Formation energy: train on Materials Project database (\( > 150,000 \) structures). Models: CGCNN (Xie & Grossman, 2018), MEGNet (Chen et al., 2019). MAE \( \approx 0.03 \) eV/atom.
Band gap: direct DFT band gap labels, corrected by scissors operator or HSE06 reference. GNN models attain MAE \( \approx 0.2 \) eV on held-out Materials Project data.
Elastic constants: predict the full \( 6 \times 6 \) stiffness tensor; requires equivariant output layers to ensure correct tensorial transformation under rotation.
Phonon properties: either train on forces across distorted configurations (enabling phonon calculation via finite differences), or predict phonon DOS directly from structure.

Example 9.4: Formation energy prediction with CGCNN.

The Crystal Graph Convolutional Neural Network (CGCNN, Xie & Grossman 2018) represents each crystal as a multigraph: nodes are atoms (feature: atomic number, electronegativity, covalent radius, …), edges are bonds (feature: distance Gaussian basis). Three convolutional layers aggregate neighbourhood information; a pooling layer averages node features; two fully-connected layers output a scalar energy. Trained on 28,046 DFT-computed structures from the Materials Project, CGCNN achieves a test MAE of 0.039 eV/atom on formation energy — dramatically better than any empirical model. This accuracy is sufficient to distinguish thermodynamically stable from metastable phases in high-throughput screening.

9.9 Connecting Classical MD, DFT, and ML

The four simulation methods covered in this course are not competing alternatives but complementary tools at different rungs of the computational hierarchy:

DFT generates the ground truth: energies, forces, electronic structure. It parametrises classical force fields and generates ML training data.
Classical MD with validated force fields reaches microsecond and micrometre scales. It provides equilibrated structures for DFT input and tests thermodynamic consistency.
ML potentials bridge the gap: DFT-accuracy potential energy surfaces at MD-accessible time scales. They are validated against both DFT (energetics) and experiment (thermodynamic and dynamical observables from MD).
kMC uses DFT-computed activation barriers to reach macroscopic time scales for rare-event kinetics: grain growth, phase transformation, thin film deposition.
MC probes thermodynamic equilibrium without dynamics: phase diagrams, adsorption isotherms, defect concentrations. It complements MD when dynamical information is not needed.

A state-of-the-art computational materials workflow might: run DFT to characterise a new alloy; fit an ML potential; run ML-MD to compute diffusion coefficients and elastic constants; use kMC with DFT barriers to compute grain boundary migration rates; and compare all results against experiment. NE 451 equips students to participate in each step of this hierarchy.

Out-of-distribution reliability. ML models are interpolators: predictions are reliable only for configurations similar to the training set. For high-strain configurations, defects not in training, or new chemistry, model predictions may fail silently. Always validate ML potential predictions against DFT for configurations critical to your conclusion. Active learning and uncertainty quantification are essential for production use in materials discovery.

Summary of Key Equations

Concept	Equation
Velocity-Verlet position	\( \mathbf{r}(t+\Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\Delta t + \tfrac{1}{2}\mathbf{a}(t)(\Delta t)^2 \)
Velocity-Verlet velocity	\( \mathbf{v}(t+\Delta t) = \mathbf{v}(t) + \tfrac{1}{2}[\mathbf{a}(t)+\mathbf{a}(t+\Delta t)]\Delta t \)
LJ potential	\( u_{\text{LJ}}(r) = 4\varepsilon[(\sigma/r)^{12} - (\sigma/r)^6] \)
Virial pressure	\( P = Nk_BT/V + \langle\sum_{i
Einstein relation	\( \text{MSD}(t) = 6Dt \)
Metropolis acceptance	\( p_{\text{acc}} = \min(1, e^{-\beta\Delta U}) \)
KS effective potential	\( v_{\text{eff}} = v_{\text{ext}} + v_H + v_{xc} \)
MLIP loss	\( \mathcal{L} = \lambda_E\sum_s(E_\theta^{(s)} - E_{\text{DFT}}^{(s)})^2/N_s^2 + \lambda_F\sum_{s,i}\\|\mathbf{F}_\theta - \mathbf{F}_{\text{DFT}}\\|^2 \)
Message passing update	\( \mathbf{h}_i^{(\ell+1)} = \phi_{\text{upd}}(\mathbf{h}_i^{(\ell)}, \sum_{j\in\mathcal{N}(i)}\phi_{\text{msg}}(\mathbf{h}_i^{(\ell)},\mathbf{h}_j^{(\ell)},\mathbf{e}_{ij})) \)
BKL time step	\( \Delta t = -\ln(u)/R, \quad R = \sum_\ell k_\ell \)