Force Fields & Molecular Dynamics¶
Classical force fields and molecular dynamics provide practical tools for large-scale atomistic simulations. With this combination of methods, you can:
Simulate material behavior up to hundreds of nanoseconds with millions of atoms
Predict thermal and mechanical properties from atomic motion
Extract engineering parameters for continuum models
The two topics are discussed together because they are often linked in practice: classical force fields define the potential energy landscape that molecular dynamics simulations explore through atomic motion, and they do it at speeds that make large-scale simulations feasible.
1. From Quantum Complexity to Classical Approximation¶
In the introduction, we learned that calculating the total energy of a material system is fundamental to understanding its behavior. The total energy determines equilibrium structures, phase stability, and mechanical properties. We also saw that the many-body problem, the quantum mechanical interaction of all electrons with each other and with nuclei, makes exact calculations impossible. We also saw that practical quantum mechanical calculations are computationally expensive for systems with many atoms. We will cover these quantum methods in more detail in a later document, but for now, let’s summarize the key challenges.
The Quantum Challenge:
One needs to solve the Schrödinger equation for all electrons simultaneously
Electrons interact with each other through Coulomb repulsion (many-body problem)
Even with DFT approximations, computational cost scales steeply with system size
DFT calculations typically limited to:
~10,000 atoms maximum
~1-10 picoseconds of dynamics (ab initio MD)
Hours, or even days, of computation time per configuration
This creates a fundamental barrier: many important engineering phenomena occur at length scales of micrometers and time scales of nanoseconds to seconds, far beyond what quantum mechanics can reach.
The Classical Solution:
For large systems and long timescales, we use classical approximations:
Force Fields (Interatomic Potentials): Replace quantum mechanical electron calculations with empirical mathematical functions
Molecular Dynamics: Simulate atomic motion using Newton’s classical equations of motion (\(F = ma\))
The strategy is conceptually simple but powerful: instead of solving for the quantum mechanical behavior of electrons, we use simpler mathematical functions that approximate how atoms interact. These functions are derived from quantum mechanics, experimental data, or both, but once established, they can be evaluated millions of times faster than solving the quantum equations.
This document explores these classical methods that enable simulation of millions of atoms over nanosecond timescales — helping to bridge the gap between quantum mechanics and continuum engineering models.
What Do We Gain?¶
The move from quantum mechanics to classical force fields provides dramatic computational advantages:
System Size Scaling: From 1,000s of atoms (DFT) → 1,000,000+ atoms (FF)
Time Scale Expansion: From picoseconds (DFT-MD) → 100s of nanoseconds (classical MD)
Accessible Engineering Phenomena:
Thermal properties: Melting transitions, thermal expansion coefficients, heat capacity, thermal conductivity
Mechanical behavior: Elastic constants, stress-strain curves, yield mechanisms, fracture toughness
Microstructure evolution: Grain growth kinetics, phase boundary migration, precipitation dynamics
Defect processes: Vacancy and interstitial diffusion, dislocation nucleation and glide, grain boundary sliding
Large-scale phenomena: Crack tip dynamics, void growth and coalescence, shock wave propagation
The computational efficiency means that parameter studies become practical: you can vary temperature, pressure, composition, or grain size and run dozens to hundreds of simulations to map out material behavior.
What Do We Lose?¶
The trade-off for computational speed is the loss of explicit quantum mechanical treatment of electrons. This imposes important limitations on what force fields and classical MD can achieve.
Electronic Structure Information
Chemical Reactions
Accuracy
Transferability
The Fundamental Trade-off¶
Classical force fields are empirical approximations to the true quantum mechanical interactions. They sacrifice completeness and universality for computational efficiency. The art of using force fields lies in:
Choosing a force field appropriate for your material and conditions
Understanding what properties the force field was designed to reproduce
Validating predictions against experiments or DFT when possible
Recognizing when you have exceeded the force field’s domain of validity
The key question that drives the rest of this module: How do we describe atomic interactions classically? What mathematical forms can capture the physics of bonding while remaining simple enough to evaluate billions of times?
2. Classical Force Fields: Mathematical Description¶
What is a Force Field?¶
A force field (also called “interatomic potential” or “potential energy function”) is a mathematical function that describes the potential energy of a system as a function of atomic positions:
Forces are derived from the potential:
The exact functional form depends on the type of material and bonding. Different classes of force fields have been developed for different materials:
2.1 Bonded Force Fields¶
For organic molecules, proteins, polymers
Most bonded force fields separate total energy into bonded and non-bonded terms:
The bonded terms describe the energy associated with chemical bonds (molecules, polymers):
The unbonded terms describe non-bonded interactions (van der Waals, electrostatics).
Examples:
AMBER (Assisted Model Building with Energy Refinement)
CHARMM (Chemistry at Harvard Macromolecular Mechanics)
OPLS (Optimized Potentials for Liquid Simulations)
GROMOS (GROningen MOlecular Simulation)
Characteristics:
Fixed bond topology (bonds do not break)
Extensive parameterization for organic chemistry
Highly accurate for biomolecules
Parameters fitted to experimental data and QM calculations
Typical Applications:
Protein folding and dynamics
Drug-receptor binding
Polymer properties
Liquid simulations
2.2 Metallic Force Fields¶
For metals and alloys
Examples:
EAM potentials for FCC metals (Al, Cu, Ni, Au)
Modified EAM (MEAM) for BCC and HCP metals
Finnis-Sinclair potentials
Embedded Atom Method (EAM):
Where:
\(F_\alpha\) = embedding function (energy to place atom in electron density)
\(\rho_\beta\) = electron density contribution from atom j
\(\phi_{\alpha\beta}\) = pair potential
Characteristics:
Captures metallic bonding (delocalized electrons)
Good for mechanical properties
Describes elastic constants, defects, surfaces
Transferable across different structures
Typical Applications:
Mechanical properties of alloys
Dislocation dynamics
Grain boundary structure
Fracture mechanics
2.3 Covalent/Semiconductor Force Fields¶
For covalent materials like Si, C, SiC
Examples:
Tersoff potential (Si, C, Ge)
Stillinger-Weber potential (Si)
REBO (Reactive Empirical Bond Order) for hydrocarbons
ReaxFF (reactive force field)
Tersoff Potential:
Where \(b_{ij}\) depends on bond angles (captures directionality of covalent bonds)
Characteristics:
Captures angular dependence of covalent bonding
Can describe bond breaking/forming (ReaxFF)
Good for semiconductors and ceramics
Typical Applications:
Silicon devices
Carbon nanostructures (nanotubes, graphene)
Chemical reactions (ReaxFF)
Ceramic materials
2.4 Ionic Force Fields¶
For ionic crystals and ceramics
Born-Mayer-Huggins Potential:
Examples:
Buckingham potential
Core-shell models (for polarizable ions)
Typical Applications:
Oxides (MgO, Al₂O₃, ZrO₂)
Solid electrolytes
Ceramic materials
Force Field Parameters¶
In addition to the functional form, force fields require a set of parameters to fill out the equations. Together, the functional form and parameters define the force field.
How are the Parameters Determined?¶
1. Fitting to Experimental Data
Elastic constants, lattice parameters
Cohesive energies
Defect formation energies
Phase transition temperatures
2. Fitting to Quantum Mechanical Calculations
DFT calculations on small systems
Energies of different configurations
Force matching to DFT forces
3. Empirical Adjustment
Trial and error refinement
Optimization to reproduce target properties
Key Consideration: Transferability
Force fields work best for systems similar to those used in parameterization
Testing on new systems is essential
Trade-off between accuracy and generality
3. Introduction to Molecular Dynamics (MD)¶
What is Molecular Dynamics?¶
Molecular Dynamics is a computational method that simulates the time evolution of a system of atoms by solving Newton’s equations of motion:
Where:
\(F_i\) = force on atom i
\(m_i\) = mass of atom i
\(r_i\) = position of atom i
\(t\) = time
The Central Role of Force Fields:
To solve Newton’s equations, we need forces on each atom. In classical MD simulations, forces come from the force field (interatomic potential) described in Section 2:
The force field \(U\) defines the physics of the system. Examples include:
EAM potential for metals (Al, Cu, Ni) → captures metallic bonding
Tersoff potential for semiconductors (Si, C) → captures covalent bonds
CHARMM/AMBER for biomolecules → captures molecular interactions
Choice of force field determines:
What materials can be simulated
Accuracy of predicted properties
Computational efficiency
Types of phenomena that can be captured
The MD Algorithm: Step by Step¶
Step 1: Initialize the System¶
Set up atomic positions (crystal structure, amorphous, etc.)
Define simulation cell and boundary conditions
Select force field appropriate for the material system
Assign initial velocities (typically drawn from a Maxwell-Boltzmann distribution at the desired temperature)
Step 2: Calculate Forces¶
Evaluate force field to compute forces on each atom
\(F_i = -\nabla_i U(r_1, r_2, ..., r_N)\)
Forces depend on positions of neighboring atoms
This is the most computationally expensive step - generally scales with the number of atoms
Cutoff distances used to limit force calculations to nearby atoms
Step 3: Integrate Equations of Motion¶
Update positions and velocities using numerical integration
Typical timestep: 0.5-5 femtoseconds (fs)
Timestep size is limited by fastest vibrations in the simulation
Step 4: Apply Constraints¶
Control temperature (thermostat) - rescale velocities or use Nosé-Hoover
Control pressure (barostat) - adjust cell volume
Apply external forces or deformations
Step 5: Repeat¶
Continue for desired simulation time
Collect data for analysis
Typical simulations: 10 ps to 100 ns
Force field evaluation dominates computational cost
Understanding MD Output: The Trajectory¶
The primary output of an MD simulation is the trajectory: a complete record of how the system evolves in time.
What is a Trajectory?¶
A trajectory stores the time evolution of the system:
Atomic positions \(\{r_i(t)\}\) for all N atoms at each timestep
Atomic velocities \(\{v_i(t)\}\) for all atoms (derived from positions or stored directly)
Simulation cell dimensions and shape (may change in NPT simulations)
Time stamps corresponding to each saved configuration
Often also: forces on atoms, energies (kinetic, potential, total), temperature, pressure
Saving this information for every timestep is not feasible in practice, and a step interval must be chosen. For a simulation with 10,000 atoms saved every 100 steps over 1 million steps, the trajectory contains 10,000 “snapshots” of the system, each recording the positions and velocities of all 10,000 atoms.
From Trajectory to Properties¶
The trajectory is raw data—you must post-process it to extract meaningful properties:
Thermodynamic properties: Average energy, pressure, volume over trajectory
Structural properties: Radial distribution functions, bond angle distributions
Mechanical properties: Stress-strain curves from deformation trajectories
Transport properties: Diffusion coefficients from mean-squared displacement
Dynamic properties: Vibrational spectra from velocity autocorrelations
The trajectory is essentially a detailed atomic-scale movie of material behavior. The time evolution of the system is recorded, an interesting events, such as a vacancy jump, a dislocation glide plane, a grain boundary migration step, may be analyzed further. This is the power of MD: not just predicting average properties, but revealing the atomic mechanisms that produce those properties.
The Trajectory is a Computational Microscope¶
The trajectory allows you to observe atomic-scale processes directly, such as watching crystals melt: You can observe how the atomic structure transitions from ordered to disordered.
4. Material Properties from MD Simulations¶
What Can We Calculate?¶
MD simulations provide access to a wealth of material properties:
Structural Properties¶
Radial distribution functions (RDF)
Coordination numbers
Bond angle distributions
Crystal structure analysis
Thermodynamic Properties¶
Temperature, pressure, volume
Heat capacity
Free energy (advanced techniques)
Phase diagrams
Mechanical Properties¶
Elastic constants (C₁₁, C₁₂, C₄₄, etc.)
Bulk modulus, shear modulus, Young’s modulus
Poisson’s ratio
Yield strength, fracture toughness
Transport Properties¶
Diffusion coefficients
Viscosity
Thermal conductivity
Ionic conductivity
Dynamic Properties¶
Vibrational density of states
Phonon dispersion
Mean squared displacement
Time correlation functions
5. Simple Application Example: Thermal Expansion of Aluminum¶
Let’s walk through a concrete example to understand how MD simulations work in practice.
The Problem¶
Goal: Calculate the coefficient of thermal expansion (CTE) for aluminum
The linear coefficient of thermal expansion is:
For volumetric expansion:
Experimental value for Al: \(\alpha_V \approx 69 \times 10^{-6}\) K-1 at 300 K
The Simulation Approach¶
Step 1: Build the Model and Select Force Field¶
Create FCC aluminum supercell (e.g., 10×10×10 unit cells with 4 atoms per unit cell = 4,000 atoms)
Apply periodic boundary conditions in all directions
Choose EAM potential for aluminum:
Why EAM? Metallic bonding requires electron density effects
Alternative: Modified EAM (MEAM)
Key consideration: Does potential reproduce experimental lattice parameter and elastic constants?
Step 2: Equilibration¶
Run NVT simulation at target temperature (e.g., 300 K)
Allow system to reach thermal equilibrium (~10-50 ps)
Monitor temperature and energy convergence
Force field quality check:
Total energy should be stable (good energy conservation in NVE test)
Temperature fluctuations should be reasonable (~10 K around target)
No structural collapse or anomalies
Step 3: NPT Production Run¶
Switch to NPT ensemble (constant pressure = 0 atm)
Run for sufficient time to get good statistics (100-500 ps)
Record volume as function of time
EAM force field enables:
Isotropic pressure control (metallic systems respond uniformly)
Cell shape fluctuations if using flexible cell
Accurate volume-temperature relationship
Step 4: Repeat at Different Temperatures¶
Perform both equilibration and production simulations at multiple temperatures (e.g., 100 K, 200 K, 300 K, 400 K, 500 K)
Each simulation gives equilibrium volume at that temperature
Temperature range limited by force field validity:
Below ~50 K: quantum effects important (FF-MD is classical)
Above ~800 K (for Al): approaching melting, ensure potential is well-parameterized
Best accuracy near temperature where potential was fitted
Step 5: Calculate CTE¶
Plot \(V(T)\) vs \(T\)
Fit linear or polynomial function
Calculate slope: \(\alpha_V = \frac{1}{V}\frac{dV}{dT}\)
Compare to experimental value: This validates the force field quality
Expected Results¶
What the output could look like¶
Temperature (K) |
Volume (ų) |
Lattice Parameter (Å) |
|---|---|---|
100 |
66,500 |
4.042 |
200 |
66,800 |
4.047 |
300 |
67,150 |
4.052 |
400 |
67,550 |
4.058 |
500 |
68,000 |
4.065 |
Analysis¶
From the table above:
Calculate \(dV/dT\) from slope of V vs T plot
At 300 K: \(\alpha_V = \frac{1}{67,150} \times \frac{67,550-66,800}{400-200} \approx 68 \times 10^{-6}\) K-1
Extending the Example¶
The general MD framework with appropriate modifications and selection of force fields can calculate:
Elastic constants: Apply strain, measure stress response (requires accurate force field)
Melting point: Heat until phase transition observed (requires force field parameterized at high T)
Diffusion coefficients: Track atomic displacements
Thermal conductivity: Apply temperature gradient, measure heat flux
Defect properties: Introduce vacancies, dislocations, grain boundaries
Return to front page: Atomic Scale Materials Modeling