Machine Learning Potentials: Bridging Accuracy and Speed

1. Introduction: The Central Dilemma¶

1.1 The Accuracy-Speed Tradeoff Revisited¶

Throughout this course, we have encountered a fundamental limitation in molecular simulation:

Quantum Mechanics (DFT):

✅ Accurate: solves the electronic Schrödinger equation
✅ Transferable: works for any chemical system
❌ Slow: scales as O(N³) with system size
❌ Limited: typically ~100 atoms, ~100 ps

Classical Force Fields (AMBER, OPLS, EAM):

✅ Fast: scales as O(N) or O(N log N)
✅ Large systems: millions of atoms, microseconds
❌ Fixed functional form: cannot capture all physics
❌ Limited transferability: requires re-parameterization for new chemistries

This creates an “accuracy gap”—there is no traditional method that is both accurate AND fast.

The Computational Ladder:
                                        
     Accuracy                          
        ↑                              
        │   CCSD(T)  ────── Gold standard (tiny systems, fs)
        │     ↑                        
        │    DFT     ────── Workhorse (small systems, ps)
        │     ↑                        
        │     ?      ────── THE GAP    
        │     ↑                        
        │  ReaxFF    ────── Reactive (medium systems, ns)
        │     ↑                        
        │    EAM     ────── Metals (large systems, ns)
        │     ↑                        
        │  Classical ────── Biomolecules (huge systems, µs)
        └──────────────────────────→ Speed

The Question: Can we fill this gap?

1.2 The Machine Learning Answer¶

Machine Learning Potentials (MLPs) represent a paradigm shift in molecular simulation. The core idea is elegantly simple:

Instead of deriving a functional form from physical principles (like LJ or harmonic bonds), we learn the potential energy surface directly from quantum mechanical data.

The Promise:

Near-DFT accuracy (errors < 1 meV/atom)
Near-classical speed (only 10-100× slower than classical MD)
1000× faster than DFT for same accuracy

The Reality Check: MLPs are powerful but not magic. They require:

Large amounts of high-quality training data (expensive DFT calculations)
Careful validation to ensure generalization
Domain expertise to avoid “garbage in, garbage out”

1.3 Historical Context: The Revolution¶

Timeline of Machine Learning Potentials:

2007 │ Behler & Parrinello publish Neural Network Potentials (NNP)
     │ → First practical MLP for condensed matter
     │
2010 │ Bartók introduces Gaussian Approximation Potential (GAP)
     │ → Kernel methods with SOAP descriptors
     │
2017 │ SchNet: First message-passing neural network for molecules
     │ → Deep learning enters the field
     │
2018 │ ANI: General-purpose NNP for organic molecules
     │ → Transferable across molecular space
     │
2021 │ NequIP: Equivariant neural networks
     │ → 1000× more data-efficient than previous methods
     │
2022 │ MACE: State-of-the-art equivariant architecture
     │ → Best accuracy with high efficiency
     │
2023 │ Foundation models (MACE-MP, CHGNet, M3GNet)
     │ → Pre-trained on massive datasets, zero-shot predictions
     │
2024 │ MLPs become mainstream in production research

We are now at a tipping point where MLPs are moving from research curiosities to standard tools.

2. The Fundamental Concepts¶

2.1 What is a Machine Learning Potential?¶

An MLP is a function that maps atomic positions to energy and forces:

E_{ML} = f_{ML}(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N; \mathbf{Z}_1, \mathbf{Z}_2, ..., \mathbf{Z}_N)

(1)

Where:

$\mathbf{r}_i$ = position of atom $i$
$\mathbf{Z}_i$ = atomic number (element type) of atom $i$
$f_{ML}$ = the learned function (neural network, kernel regression, etc.)

Forces are computed as the negative gradient:

\mathbf{F}_i = -\frac{\partial E_{ML}}{\partial \mathbf{r}_i}

(2)

Key Distinction from Classical Force Fields:

Classical: Fixed functional form with fitted parameters ( $\epsilon$ , $\sigma$ , $k_b$ , etc.)
ML: Flexible functional form learned from data (millions of parameters)

2.2 The Descriptor Problem: How to Represent Atoms¶

A neural network cannot directly take Cartesian coordinates as input. Why?

Physical Invariances Must Be Respected:

Translational Invariance: Moving the entire system by a vector $\mathbf{t}$ should not change the energy.
$E(\mathbf{r}_1 + \mathbf{t}, \mathbf{r}_2 + \mathbf{t}, ...) = E(\mathbf{r}_1, \mathbf{r}_2, ...)$
(3)
Rotational Invariance: Rotating the system by a matrix $\mathbf{R}$ should not change the energy.
$E(\mathbf{R}\mathbf{r}_1, \mathbf{R}\mathbf{r}_2, ...) = E(\mathbf{r}_1, \mathbf{r}_2, ...)$
(4)
Permutational Invariance: Swapping identical atoms should not change the energy.
$E(..., \mathbf{r}_i, ..., \mathbf{r}_j, ...) = E(..., \mathbf{r}_j, ..., \mathbf{r}_i, ...) \quad \text{if } Z_i = Z_j$
(5)

Solution: Descriptors

A descriptor is a mathematical representation of the local atomic environment that is invariant to translations, rotations, and permutations.

Raw coordinates → Descriptor → Neural Network → Energy
    (variant)      (invariant)    (learned)      (scalar)

2.3 Common Descriptor Families¶

Symmetry Functions (Behler-Parrinello)¶

The original approach from 2007 uses hand-crafted functions:

Radial Symmetry Functions (G2):

G_i^{(2)} = \sum_{j \neq i} e^{-\eta(r_{ij} - R_s)^2} \cdot f_c(r_{ij})

(6)

This measures “how many neighbors are at distance $R_s$ from atom $i$ ?”

Angular Symmetry Functions (G4):

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

(7)

This measures “what is the angular distribution of neighbors around atom $i$ ?”

Cutoff Function (smoothly goes to zero):

f_c(r) = \begin{cases} 0.5 \left[1 + \cos\left(\frac{\pi r}{r_c}\right)\right] & r \leq r_c \\ 0 & r > r_c \end{cases}

(8)

Total Descriptor Vector: For each atom, we compute many symmetry functions with different parameters ( $\eta$ , $R_s$ , $\zeta$ , $\lambda$ ), creating a vector of ~50-200 values that encodes the local environment.

SOAP (Smooth Overlap of Atomic Positions)¶

A more principled approach based on atomic density:

Represent neighbors as Gaussian densities centered on each atom
Expand this density in spherical harmonics and radial basis functions
Take the power spectrum (rotationally invariant)

p_{nn'l}^i = \pi \sqrt{\frac{8}{2l+1}} \sum_m c_{nlm}^{i*} c_{n'lm}^i

(9)

Advantages:

Smooth and differentiable everywhere
Systematic: can increase resolution by adding more basis functions
Well-suited for kernel methods (GAP)

Graph Representations (Modern Approaches)¶

Modern MLPs represent atoms as nodes in a graph:

Nodes: Atoms (features = element embedding)
Edges: Interatomic distances and directions

Message-passing neural networks update node features by aggregating information from neighbors:

\mathbf{h}_i^{(l+1)} = \phi\left(\mathbf{h}_i^{(l)}, \sum_{j \in \mathcal{N}(i)} \psi(\mathbf{h}_i^{(l)}, \mathbf{h}_j^{(l)}, \mathbf{e}_{ij})\right)

(10)

This learns descriptors rather than hand-crafting them.

2.4 Interactive Example: Visualizing Symmetry Functions¶

import numpy as np
import matplotlib.pyplot as plt

def cutoff_function(r, r_c):
    """Cosine cutoff function"""
    result = np.zeros_like(r)
    mask = r <= r_c
    result[mask] = 0.5 * (1 + np.cos(np.pi * r[mask] / r_c))
    return result

def G2_symmetry_function(r, eta, R_s, r_c):
    """Radial symmetry function G2"""
    return np.exp(-eta * (r - R_s)**2) * cutoff_function(r, r_c)

# Parameters
r_c = 6.0  # Cutoff radius (Angstroms)
r = np.linspace(0.5, 7.0, 500)

# Different G2 functions with varying parameters
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# --- Plot 1: Effect of R_s (peak position) ---
ax1 = axes[0, 0]
eta_fixed = 0.5
R_s_values = [1.5, 2.5, 3.5, 4.5]
colors = plt.cm.viridis(np.linspace(0, 0.8, len(R_s_values)))

for R_s, color in zip(R_s_values, colors):
    G2 = G2_symmetry_function(r, eta_fixed, R_s, r_c)
    ax1.plot(r, G2, color=color, linewidth=2.5, label=f'R_s = {R_s} Å')

ax1.axvline(r_c, color='red', linestyle='--', alpha=0.5, label=f'Cutoff = {r_c} Å')
ax1.set_xlabel('Distance r (Å)', fontsize=12)
ax1.set_ylabel('G2 Value', fontsize=12)
ax1.set_title(f'Radial Symmetry Functions (η = {eta_fixed})', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)
ax1.set_xlim(0.5, 7)

# --- Plot 2: Effect of eta (width) ---
ax2 = axes[0, 1]
R_s_fixed = 2.5
eta_values = [0.1, 0.5, 1.0, 2.0]
colors = plt.cm.plasma(np.linspace(0, 0.8, len(eta_values)))

for eta, color in zip(eta_values, colors):
    G2 = G2_symmetry_function(r, eta, R_s_fixed, r_c)
    ax2.plot(r, G2, color=color, linewidth=2.5, label=f'η = {eta}')

ax2.axvline(r_c, color='red', linestyle='--', alpha=0.5, label=f'Cutoff = {r_c} Å')
ax2.set_xlabel('Distance r (Å)', fontsize=12)
ax2.set_ylabel('G2 Value', fontsize=12)
ax2.set_title(f'Effect of Width Parameter η (R_s = {R_s_fixed} Å)', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)
ax2.set_xlim(0.5, 7)

# --- Plot 3: Cutoff function ---
ax3 = axes[1, 0]
cutoff_values = [4.0, 5.0, 6.0, 7.0]
colors = plt.cm.coolwarm(np.linspace(0.2, 0.8, len(cutoff_values)))

for rc, color in zip(cutoff_values, colors):
    fc = cutoff_function(r, rc)
    ax3.plot(r, fc, color=color, linewidth=2.5, label=f'r_c = {rc} Å')

ax3.set_xlabel('Distance r (Å)', fontsize=12)
ax3.set_ylabel('f_c(r)', fontsize=12)
ax3.set_title('Cosine Cutoff Function', fontsize=13, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(alpha=0.3)
ax3.set_xlim(0.5, 7)

# --- Plot 4: Complete descriptor set ---
ax4 = axes[1, 1]

# A practical descriptor set
descriptor_params = [
    (0.3, 1.5), (0.3, 2.5), (0.3, 3.5), (0.3, 4.5),
    (1.0, 2.0), (1.0, 3.0), (1.0, 4.0),
]

colors = plt.cm.rainbow(np.linspace(0, 1, len(descriptor_params)))
for (eta, R_s), color in zip(descriptor_params, colors):
    G2 = G2_symmetry_function(r, eta, R_s, r_c)
    ax4.plot(r, G2, color=color, linewidth=2, alpha=0.8, 
             label=f'η={eta}, Rs={R_s}')

ax4.axvline(r_c, color='red', linestyle='--', alpha=0.5)
ax4.set_xlabel('Distance r (Å)', fontsize=12)
ax4.set_ylabel('G2 Value', fontsize=12)
ax4.set_title('Complete Descriptor Set (Multiple G2)', fontsize=13, fontweight='bold')
ax4.legend(fontsize=8, ncol=2)
ax4.grid(alpha=0.3)
ax4.set_xlim(0.5, 7)

plt.tight_layout()
plt.savefig('symmetry_functions.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("SYMMETRY FUNCTION INTERPRETATION:")
print("="*60)
print("Each G2 function is a 'sensor' for a specific distance range:")
print(f"  • Small η (wide) → Counts all neighbors broadly")
print(f"  • Large η (narrow) → Sensitive to specific distance")
print(f"  • R_s → Peak sensitivity location")
print(f"\nThe complete set of G2 values forms a 'fingerprint' of")
print(f"the local environment that is:")
print(f"  ✓ Translationally invariant (depends only on distances)")
print(f"  ✓ Rotationally invariant (depends only on |r_ij|)")
print(f"  ✓ Permutationally invariant (sum over all neighbors)")
print("="*60)

3. Major MLP Architectures¶

3.1 Behler-Parrinello Neural Network Potentials (2007)¶

The architecture that started it all.

Key Idea: Atomic Decomposition

The total energy is a sum of atomic contributions:

E_{total} = \sum_{i=1}^{N} E_i(\mathbf{G}_i)

(11)

Where $\mathbf{G}_i$ is the vector of symmetry functions for atom $i$ , and $E_i$ is computed by a neural network.

Architecture:

For each atom i:
    
    Local environment → Symmetry Functions → Neural Network → Atomic Energy
         (r_ij)              (G2, G4)           (MLP)           E_i
                                                                 │
                                                                 ↓
                                                    Sum over atoms: E_total

Neural Network Structure:

Typically 2-3 hidden layers
20-50 neurons per layer
tanh or softplus activation
Separate networks for each element type

Strengths:

Well-understood, extensively validated
Efficient (linear scaling with system size)
Good for homogeneous systems (bulk metals, simple alloys)

Limitations:

Requires hand-crafted descriptors (choosing η, R_s, etc.)
Angular functions expensive for many neighbors
Limited to ~3-4 element types efficiently

3.2 Gaussian Approximation Potential (GAP)¶

A kernel-based approach using SOAP descriptors.

Key Idea: Kernel Regression

Instead of training a neural network, GAP uses kernel ridge regression:

E_i = \sum_{k=1}^{M} \alpha_k \kappa(\mathbf{p}_i, \mathbf{p}_k)

(12)

Where:

$\mathbf{p}_i$ = SOAP descriptor of atom $i$
$\mathbf{p}_k$ = descriptors of training configurations
$\kappa$ = kernel function (measures similarity)
$\alpha_k$ = learned coefficients

The SOAP Kernel:

\kappa(\mathbf{p}_i, \mathbf{p}_j) = \left( \frac{\mathbf{p}_i \cdot \mathbf{p}_j}{|\mathbf{p}_i| |\mathbf{p}_j|} \right)^\zeta

(13)

Strengths:

Excellent accuracy with limited training data
Uncertainty quantification (knows when it doesn’t know)
Mature, well-validated software (QUIP/GAP)

Limitations:

Scales poorly: evaluation cost O(M) where M = number of sparse points
Memory-intensive for large training sets
Requires careful sparsification for practical use

3.3 Message-Passing Neural Networks (SchNet, DimeNet)¶

Modern deep learning approaches that learn descriptors.

Key Idea: Learn Everything End-to-End

Instead of hand-crafted descriptors, use graph neural networks to learn atomic representations:

Architecture:

Initial embedding:        h_i^(0) = embed(Z_i)    (element-specific vector)
                               ↓
Message passing (L layers):  
    For l = 1 to L:
        m_ij = Message(h_i^(l-1), h_j^(l-1), r_ij)  (information from neighbor j)
        h_i^(l) = Update(h_i^(l-1), Σ_j m_ij)       (aggregate messages)
                               ↓
Readout:                  E_i = Readout(h_i^(L))   (atomic energy)
                               ↓
Sum:                      E_total = Σ_i E_i

SchNet (2017):

Uses continuous-filter convolutional layers
Filters depend on interatomic distances
First to show transferability across molecules

DimeNet (2020):

Incorporates angles explicitly (2-body and 3-body)
Higher accuracy than SchNet
More expensive computationally

PaiNN (2021):

Equivariant message passing
Both scalar and vector features
Good balance of accuracy and speed

3.4 Equivariant Neural Networks (NequIP, MACE)¶

The current state-of-the-art.

Key Insight: Equivariance, not Invariance

Previous methods used invariant descriptors (scalars). But forces are vectors!

Invariance: Output unchanged under transformation

f(\mathbf{R}\mathbf{x}) = f(\mathbf{x})

(14)

Equivariance: Output transforms predictably

f(\mathbf{R}\mathbf{x}) = \mathbf{R}f(\mathbf{x})

(15)

If we learn equivariant features, we get forces “for free” by taking gradients.

NequIP (2021):

Uses E(3)-equivariant neural networks
Features are spherical tensors, not just scalars
100-1000× more data-efficient than invariant networks

MACE (2022):

Many-body equivariant message passing
Higher body-order correlations without explicit enumeration
Current best accuracy on standard benchmarks
Foundation models available (MACE-MP-0)

# Conceptual comparison of architectures

import matplotlib.pyplot as plt
import numpy as np

# Comparison data (approximate, based on literature)
architectures = ['BP-NNP\n(2007)', 'GAP/SOAP\n(2010)', 'SchNet\n(2017)', 
                 'NequIP\n(2021)', 'MACE\n(2022)']

# Error on MD17 benchmark (meV/atom for ethanol)
mae_error = [2.0, 1.5, 1.0, 0.3, 0.15]

# Training data needed for 1 meV/atom (relative, log scale)
data_needed = [10000, 5000, 2000, 200, 100]

# Inference speed (relative to classical, lower is faster)
speed = [50, 200, 100, 80, 60]

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# --- Plot 1: Accuracy ---
ax1 = axes[0]
colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(architectures)))
bars1 = ax1.bar(architectures, mae_error, color=colors, edgecolor='black', linewidth=1.5)
ax1.set_ylabel('MAE (meV/atom)', fontsize=11)
ax1.set_title('Accuracy (Lower is Better)', fontsize=12, fontweight='bold')
ax1.set_ylim(0, 2.5)
for bar, val in zip(bars1, mae_error):
    ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, 
             f'{val}', ha='center', fontsize=10)

# --- Plot 2: Data Efficiency ---
ax2 = axes[1]
colors = plt.cm.Greens(np.linspace(0.4, 0.9, len(architectures)))
bars2 = ax2.bar(architectures, data_needed, color=colors, edgecolor='black', linewidth=1.5)
ax2.set_ylabel('Training Configs Needed', fontsize=11)
ax2.set_title('Data Efficiency (Lower is Better)', fontsize=12, fontweight='bold')
ax2.set_yscale('log')
for bar, val in zip(bars2, data_needed):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() * 1.2,
             f'{val}', ha='center', fontsize=10)

# --- Plot 3: Speed ---
ax3 = axes[2]
colors = plt.cm.Oranges(np.linspace(0.4, 0.9, len(architectures)))
bars3 = ax3.bar(architectures, speed, color=colors, edgecolor='black', linewidth=1.5)
ax3.set_ylabel('Slowdown vs Classical MD', fontsize=11)
ax3.set_title('Speed (Lower is Better)', fontsize=12, fontweight='bold')
for bar, val in zip(bars3, speed):
    ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5,
             f'{val}×', ha='center', fontsize=10)

plt.tight_layout()
plt.savefig('architecture_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("ARCHITECTURE EVOLUTION:")
print("="*60)
print("2007 BP-NNP → First practical MLP, hand-crafted descriptors")
print("2010 GAP    → Kernel methods, better accuracy, slower")
print("2017 SchNet → Deep learning, learns descriptors automatically")
print("2021 NequIP → Equivariance breakthrough, 100× less data needed")
print("2022 MACE   → Many-body equivariance, current state-of-the-art")
print("="*60)

4. Training Machine Learning Potentials¶

4.1 The Training Workflow¶

Training an MLP requires a complete pipeline:

1. Generate Training Data
   ├── Run DFT calculations on diverse structures
   ├── Extract energies, forces, (stresses)
   └── Store in standard format (XYZ, ASE database)
           ↓
2. Prepare Dataset
   ├── Split: training / validation / test
   ├── Compute statistics (mean, std for normalization)
   └── Convert to ML framework format (PyTorch, TensorFlow)
           ↓
3. Train Model
   ├── Define architecture (layers, features, cutoff)
   ├── Set loss function (energy + force weights)
   ├── Optimize with gradient descent (Adam, SGD)
   └── Monitor validation loss (early stopping)
           ↓
4. Validate Model
   ├── Evaluate on held-out test set
   ├── Run MD and check physical properties
   └── Compare to experimental/reference data
           ↓
5. Deploy
   ├── Export model to LAMMPS/ASE format
   └── Run production simulations

4.2 Training Data: Quality Over Quantity¶

Where Does Training Data Come From?

Most commonly: Density Functional Theory (DFT) calculations

What Makes Good Training Data?

Diversity: Cover the configuration space you want to explore
- Bulk structures at different densities
- Surfaces, interfaces, defects
- High-temperature / high-pressure configurations
- Reaction pathways (if applicable)
Relevance: Include structures similar to your target simulations
- Don’t train on water if you want to simulate metals
- Include relevant temperature ranges
Accuracy: Use consistent, high-quality reference calculations
- Same DFT functional throughout
- Converged basis sets and k-points
- Check for DFT artifacts (self-interaction, dispersion)

Data Generation Strategies:

# Example: Generating diverse training structures

import numpy as np

# Strategy 1: Rattled structures (add random displacements)
def rattle_structure(positions, amplitude=0.1):
    """Add random noise to atomic positions"""
    noise = np.random.randn(*positions.shape) * amplitude
    return positions + noise

# Strategy 2: Strained structures (scale cell)
def apply_strain(cell, strain_matrix):
    """Apply strain tensor to unit cell"""
    return np.dot(cell, np.eye(3) + strain_matrix)

# Strategy 3: MD snapshots (sample thermal fluctuations)
# Run short MD with approximate potential, extract snapshots

# Strategy 4: Active learning (iterative refinement)
# 1. Train initial model on small dataset
# 2. Run MD with current model
# 3. Identify high-uncertainty configurations
# 4. Calculate DFT for those configurations
# 5. Add to training set, retrain
# 6. Repeat until converged

4.3 The Loss Function¶

MLPs are trained to minimize the difference between predicted and reference values:

\mathcal{L} = w_E \sum_{\text{configs}} (E_{pred} - E_{ref})^2 + w_F \sum_{\text{atoms}} |\mathbf{F}_{pred} - \mathbf{F}_{ref}|^2 + w_S \sum_{\text{configs}} |\sigma_{pred} - \sigma_{ref}|^2

(16)

Where:

$w_E, w_F, w_S$ = weights for energy, forces, stress
Typically $w_F \gg w_E$ because forces are more informative

Why Include Forces?

Forces provide $3N$ data points per configuration vs. 1 for energy:

100 atoms → 300 force components, but only 1 energy
Forces constrain the potential energy surface more tightly

Typical Weight Ratios:

Energy: 1.0
Forces: 100.0 - 1000.0 (per atom)
Stress: 0.1 - 1.0 (optional)

4.4 Avoiding Overfitting¶

The Danger: MLPs have millions of parameters and can memorize training data.

Solutions:

Train/Validation/Test Split:
- 80% training / 10% validation / 10% test
- Monitor validation loss during training
- Stop when validation loss stops improving (early stopping)
Regularization:
- Weight decay (L2 penalty)
- Dropout (less common for MLPs)
- Noise injection on inputs
Data Augmentation:
- Random rotations (if not using equivariant network)
- Atom permutations
Ensemble Methods:
- Train multiple models with different random seeds
- Average predictions (reduces variance)
- Disagreement = uncertainty estimate

# Visualization: Learning curves

import matplotlib.pyplot as plt
import numpy as np

# Mock training curves
epochs = np.arange(1, 201)

# Good training (converged, no overfitting)
train_loss_good = 10 * np.exp(-epochs / 30) + 0.5 + 0.1 * np.random.randn(len(epochs))
val_loss_good = 10 * np.exp(-epochs / 35) + 0.6 + 0.15 * np.random.randn(len(epochs))

# Overfitting (validation diverges)
train_loss_overfit = 10 * np.exp(-epochs / 20) + 0.1 * np.exp(-epochs / 100) + 0.05 * np.random.randn(len(epochs))
val_loss_overfit = 10 * np.exp(-epochs / 30) + 0.3 + 0.02 * epochs + 0.2 * np.random.randn(len(epochs))
val_loss_overfit = np.maximum(val_loss_overfit, 0.5)  # Keep positive

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# --- Good training ---
ax1.plot(epochs, train_loss_good, 'b-', linewidth=2, label='Training Loss', alpha=0.8)
ax1.plot(epochs, val_loss_good, 'r-', linewidth=2, label='Validation Loss', alpha=0.8)
ax1.axhline(y=0.6, color='green', linestyle='--', alpha=0.5, label='Converged')
ax1.set_xlabel('Epoch', fontsize=12)
ax1.set_ylabel('Loss', fontsize=12)
ax1.set_title('Good Training: Converged', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.set_yscale('log')
ax1.set_ylim(0.3, 15)
ax1.grid(alpha=0.3)

# --- Overfitting ---
ax2.plot(epochs, train_loss_overfit, 'b-', linewidth=2, label='Training Loss', alpha=0.8)
ax2.plot(epochs, val_loss_overfit, 'r-', linewidth=2, label='Validation Loss', alpha=0.8)
ax2.axvline(x=60, color='orange', linestyle='--', linewidth=2, label='Stop Here!')
ax2.annotate('Overfitting starts', xy=(60, 3), xytext=(100, 6),
             fontsize=11, arrowprops=dict(arrowstyle='->', color='orange', lw=2))
ax2.set_xlabel('Epoch', fontsize=12)
ax2.set_ylabel('Loss', fontsize=12)
ax2.set_title('Overfitting: Validation Diverges', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.set_yscale('log')
ax2.set_ylim(0.1, 15)
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('learning_curves.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("TRAINING BEST PRACTICES:")
print("="*60)
print("1. Always monitor validation loss, not just training loss")
print("2. Use early stopping: stop when validation loss plateaus")
print("3. Save checkpoints frequently")
print("4. Compare multiple training runs (different seeds)")
print("5. Test on truly held-out data before production")
print("="*60)

4.5 Error Metrics¶

Standard metrics for evaluating MLP quality:

Energy Errors:

MAE (Mean Absolute Error): $\frac{1}{N}\sum |E_{pred} - E_{ref}|$
RMSE (Root Mean Square Error): $\sqrt{\frac{1}{N}\sum (E_{pred} - E_{ref})^2}$
Usually reported per atom: “meV/atom”

Force Errors:

MAE on force components: $\frac{1}{3N_{atoms}}\sum |F_{pred} - F_{ref}|$
Usually reported in “meV/Å” or “eV/Å”

Benchmark Targets:

Property	Chemical Accuracy	Good MLP	Acceptable
Energy	1 meV/atom	<5 meV/atom	<20 meV/atom
Forces	10 meV/Å	<50 meV/Å	<100 meV/Å

5. Active Learning: Smart Data Generation¶

5.1 The Problem with Random Sampling¶

Naively generating training data is wasteful:

Most of phase space is irrelevant to your simulation
You’ll spend DFT compute on easy configurations
Important rare events may be missed

Active Learning solves this by letting the model guide data generation.

5.2 The Active Learning Loop¶

                    ┌─────────────────────────┐
                    │                         │
                    ▼                         │
    ┌─────────────────────────┐               │
    │ 1. Train MLP on         │               │
    │    current dataset      │               │
    └───────────┬─────────────┘               │
                │                             │
                ▼                             │
    ┌─────────────────────────┐               │
    │ 2. Run MD with          │               │
    │    current MLP          │               │
    └───────────┬─────────────┘               │
                │                             │
                ▼                             │
    ┌─────────────────────────┐               │
    │ 3. Identify uncertain   │               │
    │    configurations       │               │
    └───────────┬─────────────┘               │
                │                             │
                ▼                             │
    ┌─────────────────────────┐               │
    │ 4. Run DFT on           │               │
    │    selected configs     │               │
    └───────────┬─────────────┘               │
                │                             │
                ▼                             │
    ┌─────────────────────────┐               │
    │ 5. Add to training set  │───────────────┘
    └─────────────────────────┘
    
    Repeat until convergence

5.3 Uncertainty Quantification¶

How do we know when the model is uncertain?

Method 1: Committee Disagreement

Train an ensemble of models (different random seeds)
Measure variance in predictions
High variance = high uncertainty

\sigma_E = \sqrt{\frac{1}{M}\sum_{m=1}^M (E_m - \bar{E})^2}

(17)

Method 2: Query by Committee

If committee members disagree, the model is extrapolating
Select configurations with highest disagreement

Method 3: Bayesian Methods (GAP)

GAP naturally provides uncertainty estimates
Epistemic uncertainty from posterior distribution

# Illustration of committee disagreement

import numpy as np
import matplotlib.pyplot as plt

# Mock 1D potential energy surface
x = np.linspace(-3, 3, 200)
true_energy = np.sin(2*x) + 0.3 * x**2 - 2

# Training data (sparse)
x_train = np.array([-2.5, -1.5, -0.5, 0.5, 1.5])
y_train = np.sin(2*x_train) + 0.3 * x_train**2 - 2 + 0.1*np.random.randn(len(x_train))

# Committee predictions (5 models with different fits)
np.random.seed(42)
committee_predictions = []
for i in range(5):
    # Each model fits slightly differently
    noise = 0.2 * np.sin(3*x + i) + 0.1 * np.random.randn(len(x))
    pred = true_energy + noise
    # Extrapolation region gets worse
    pred[x < -2] += 0.3 * (i - 2) * (x[x < -2] + 2)**2
    pred[x > 2] += 0.3 * (i - 2) * (x[x > 2] - 2)**2
    committee_predictions.append(pred)

committee_predictions = np.array(committee_predictions)
mean_pred = np.mean(committee_predictions, axis=0)
std_pred = np.std(committee_predictions, axis=0)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# --- Plot 1: Committee predictions ---
for i, pred in enumerate(committee_predictions):
    ax1.plot(x, pred, alpha=0.5, linewidth=1.5, label=f'Model {i+1}')
ax1.plot(x, true_energy, 'k--', linewidth=2.5, label='True PES')
ax1.scatter(x_train, y_train, s=100, c='red', zorder=5, edgecolors='black', 
            linewidths=2, label='Training Data')
ax1.set_xlabel('Configuration x', fontsize=12)
ax1.set_ylabel('Energy', fontsize=12)
ax1.set_title('Committee of 5 Models', fontsize=13, fontweight='bold')
ax1.legend(fontsize=9, ncol=2)
ax1.set_xlim(-3, 3)
ax1.grid(alpha=0.3)

# --- Plot 2: Uncertainty ---
ax2.fill_between(x, mean_pred - 2*std_pred, mean_pred + 2*std_pred, 
                  alpha=0.3, color='blue', label='±2σ uncertainty')
ax2.plot(x, mean_pred, 'b-', linewidth=2.5, label='Mean prediction')
ax2.plot(x, true_energy, 'k--', linewidth=2.5, label='True PES')
ax2.scatter(x_train, y_train, s=100, c='red', zorder=5, edgecolors='black', 
            linewidths=2, label='Training Data')

# Highlight high-uncertainty regions
high_uncertainty = std_pred > 0.15
ax2.fill_between(x, ax2.get_ylim()[0], ax2.get_ylim()[1], 
                  where=high_uncertainty, alpha=0.1, color='red')

ax2.set_xlabel('Configuration x', fontsize=12)
ax2.set_ylabel('Energy', fontsize=12)
ax2.set_title('Uncertainty-Guided Selection', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.set_xlim(-3, 3)
ax2.set_ylim(-4, 2)
ax2.grid(alpha=0.3)

# Annotate
ax2.annotate('Add training\ndata here!', xy=(-2.7, 0.5), xytext=(-2.2, 1.5),
             fontsize=11, fontweight='bold', color='red',
             arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax2.annotate('Add training\ndata here!', xy=(2.7, 0.5), xytext=(2.2, 1.5),
             fontsize=11, fontweight='bold', color='red',
             arrowprops=dict(arrowstyle='->', color='red', lw=2))

plt.tight_layout()
plt.savefig('active_learning.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("ACTIVE LEARNING WORKFLOW:")
print("="*60)
print("1. Start with small, diverse training set (~100-1000 configs)")
print("2. Train initial model")
print("3. Run MD simulation with trained model")
print("4. Monitor committee disagreement during MD")
print("5. When disagreement exceeds threshold:")
print("   → Pause MD")
print("   → Run DFT on that configuration")
print("   → Add to training set")
print("   → Optionally retrain model")
print("6. Continue MD")
print("7. Repeat until no new high-uncertainty configs appear")
print("="*60)

6. Using MLPs in LAMMPS¶

6.1 Available Interfaces¶

LAMMPS supports multiple MLP frameworks:

Package	Architecture	Command	Notes
ML-PACE	ACE (Atomic Cluster Expansion)	`pair_style pace`	Fast, flexible
ML-SNAP	SNAP (bispectrum)	`pair_style snap`	Good for metals
n2p2	Behler-Parrinello NNP	`pair_style nnp`	Well-documented
DeePMD-kit	Deep Potential	`pair_style deepmd`	GPU-accelerated
MACE	Equivariant GNN	`pair_style mace`	State-of-the-art

6.2 Example: Running with a Pre-trained Model¶

For foundation models like MACE-MP-0, you can download pre-trained weights:

# LAMMPS input file: mace_example.in

# --- Initialization ---
units           metal
atom_style      atomic
boundary        p p p

# --- Read structure ---
read_data       silicon_bulk.data

# --- Define MLP potential ---
pair_style      mace no_domain_decomposition
pair_coeff      * * mace_mp_0.model Si

# --- Settings ---
neighbor        2.0 bin
neigh_modify    every 1 delay 0 check yes

# --- Thermodynamic output ---
thermo          100
thermo_style    custom step temp pe ke etotal press vol

# --- Equilibration (NVT) ---
velocity        all create 300.0 12345 mom yes rot yes
fix             1 all nvt temp 300.0 300.0 0.1
timestep        0.001   # 1 fs

run             10000

# --- Production (NVT) ---
dump            1 all custom 100 trajectory.dump id type x y z vx vy vz
run             50000

6.3 Example: Training Your Own Model with MACE¶

# Python script: train_mace.py
# This is a simplified example - see MACE documentation for full options

from mace.cli.run_train import main as mace_train
import subprocess

# Prepare training configuration
config = """
# MACE training configuration
name: my_water_model
seed: 42

# Model architecture
model: MACE
hidden_irreps: 128x0e + 128x1o
r_max: 5.0
num_radial_basis: 8
num_cutoff_basis: 5
max_ell: 3
correlation: 3
num_interactions: 2

# Training data
train_file: water_train.xyz
valid_file: water_valid.xyz
energy_key: REF_energy
forces_key: REF_forces

# Training parameters
batch_size: 10
max_num_epochs: 500
learning_rate: 0.01
scheduler: ReduceLROnPlateau
patience: 50

# Loss weights
energy_weight: 1.0
forces_weight: 100.0

# Output
model_dir: ./trained_models
"""

# Write config file
with open('mace_config.yaml', 'w') as f:
    f.write(config)

# Train model (command line)
print("To train, run:")
print("mace_run_train --config mace_config.yaml")

6.4 Computational Considerations¶

GPU Acceleration: Most modern MLPs benefit enormously from GPUs:

DeePMD-kit: 10-50× speedup on GPU
MACE: 5-20× speedup on GPU
CPU-only is viable for small systems

Memory Requirements:

Model weights: typically 10-100 MB
Per-atom buffers: depends on cutoff and neighbor count
Batch size affects GPU memory

Scaling:

MLPs scale as O(N) with system size (like classical potentials)
But with larger prefactor (typically 10-100× slower per step)
Parallelization over atoms is efficient

# Performance comparison visualization

import matplotlib.pyplot as plt
import numpy as np

# System sizes
N_atoms = np.array([100, 500, 1000, 5000, 10000, 50000])

# Time per step (relative, arbitrary units)
time_classical = 0.001 * N_atoms  # Linear scaling
time_dft = 0.1 * N_atoms**3  # Cubic scaling (approximate)
time_mlp = 0.05 * N_atoms   # Linear but higher prefactor

# Cap DFT at practical limit
time_dft[N_atoms > 500] = np.nan

fig, ax = plt.subplots(figsize=(10, 6))

ax.loglog(N_atoms, time_classical, 'g-o', linewidth=2.5, markersize=8, 
          label='Classical FF (OPLS)')
ax.loglog(N_atoms, time_mlp, 'b-s', linewidth=2.5, markersize=8,
          label='MLP (MACE)')
ax.loglog(N_atoms[~np.isnan(time_dft)], time_dft[~np.isnan(time_dft)], 
          'r-^', linewidth=2.5, markersize=8, label='DFT (AIMD)')

# Annotate regions
ax.fill_between([100, 500], [1e-3, 1e-3], [1e4, 1e4], alpha=0.1, color='red')
ax.text(200, 1e2, 'DFT\nFeasible', fontsize=11, ha='center', color='red')

ax.fill_between([1000, 50000], [1e-3, 1e-3], [1e4, 1e4], alpha=0.1, color='blue')
ax.text(5000, 1e2, 'MLP\nSweet Spot', fontsize=11, ha='center', color='blue')

ax.set_xlabel('Number of Atoms', fontsize=12)
ax.set_ylabel('Time per MD Step (arb. units)', fontsize=12)
ax.set_title('Computational Cost: Classical vs MLP vs DFT', fontsize=13, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(alpha=0.3, which='both')
ax.set_xlim(80, 70000)
ax.set_ylim(5e-2, 1e5)

plt.tight_layout()
plt.savefig('cost_comparison.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("COMPUTATIONAL COST SUMMARY:")
print("="*60)
print("System Size   | Classical  | MLP       | DFT")
print("--------------|------------|-----------|----------")
print("100 atoms     | 0.001s     | 0.1s      | 60s")
print("1,000 atoms   | 0.01s      | 1s        | Impractical")
print("10,000 atoms  | 0.1s       | 10s       | Impossible")
print("="*60)
print("\nMLP fills the gap: DFT accuracy at near-classical speed")
print("for systems of 1,000-100,000 atoms.")
print("="*60)

7. Limitations and Best Practices¶

7.1 The Extrapolation Problem¶

The Fundamental Limitation: MLPs are interpolators, not extrapolators. They work well inside the training distribution but can fail catastrophically outside it.

Warning Signs:

Unphysical forces (atoms flying apart)
Energy drifting during MD
Melting point far from experiment
Negative phonon frequencies

Solutions:

Comprehensive training data: Cover relevant phase space
Active learning: Catch extrapolation during simulation
Uncertainty monitoring: Flag suspicious configurations
Validation: Always compare to experiment/DFT

7.2 Common Pitfalls¶

Pitfall	Symptom	Solution
Training data not diverse enough	Model fails at new T, P	Add more configurations
Overfitting	Low train error, high test error	Regularization, more data
Cutoff too short	Missing long-range interactions	Increase cutoff, add corrections
Wrong reference method	Systematic bias	Use consistent DFT settings
Ignoring periodicity	Surface artifacts in bulk	Check boundary conditions

7.3 Validation Checklist¶

Before trusting MLP results for production:

Static Properties:

Lattice constants within 1% of experiment
Elastic constants within 10%
Phonon dispersion qualitatively correct
Surface energies reasonable

Dynamic Properties:

Melting point within 50-100 K
Diffusion coefficients correct order of magnitude
Radial distribution functions match experiment
No energy drift in long NVE simulations

Stability:

MD runs for >1 ns without crashes
No atoms with unphysical velocities
Temperature/pressure fluctuations reasonable

8. Foundation Models: The New Frontier¶

8.1 What are Foundation Models?¶

Traditional MLPs: Train a new model for each system (water model, silicon model, etc.)

Foundation Models: Train once on massive, diverse datasets, then use for many systems without retraining.

Key Examples (2023-2024):

MACE-MP-0: Trained on Materials Project data (~150,000 materials)
CHGNet: Trained on Materials Project with charge information
M3GNet: Similar scope, different architecture
GNoME: Google’s model for materials discovery

8.2 When to Use Foundation Models¶

Good Use Cases:

✅ Quick screening of many materials
✅ Initial exploration before expensive DFT
✅ Systems similar to training data (inorganic crystals)
✅ Rapid prototyping of simulations

Proceed with Caution:

⚠️ Organic molecules (less training data)
⚠️ Extreme conditions (high T, high P)
⚠️ Defects and interfaces (underrepresented)
⚠️ Quantitative predictions (validate first!)

Not Recommended:

❌ Systems far from training distribution
❌ When high accuracy is critical
❌ Novel chemistries not in training data

8.3 Fine-Tuning¶

Foundation models can be adapted to your specific system:

Foundation Model (general)
         │
         │ Fine-tune on 100-1000 system-specific DFT calculations
         ▼
Specialized Model (your system, higher accuracy)

This combines the broad knowledge of the foundation model with system-specific accuracy.

9. Practical Workflow Summary¶

9.1 Decision Tree: Which Potential to Use?¶

Start
  │
  ├── Is my system well-described by classical FF?
  │     │
  │     ├── Yes → Use AMBER/CHARMM/OPLS (fast, validated)
  │     │
  │     └── No → Continue
  │
  ├── Do I need reactive chemistry?
  │     │
  │     ├── Yes → Consider ReaxFF or MLP
  │     │
  │     └── No → Continue
  │
  ├── Is my system a simple metal/alloy?
  │     │
  │     ├── Yes → Use EAM (fast, accurate for metals)
  │     │
  │     └── No → Continue
  │
  ├── Is my system in the Materials Project database?
  │     │
  │     ├── Yes → Try MACE-MP-0 / CHGNet (zero-shot)
  │     │         Validate before production!
  │     │
  │     └── No → Continue
  │
  ├── Do I have DFT data for my system?
  │     │
  │     ├── Yes → Train custom MLP (NequIP, MACE)
  │     │
  │     └── No → Generate DFT data first, then train MLP
  │
  └── Consider: Is the MLP effort worth it?
        - Small system, short time → Just use DFT
        - Large system, long time → MLP is justified

9.2 Recommended Workflow¶

Start Simple:
- Check if classical FF exists for your system
- Try foundation model (zero-shot)
- Validate against known properties
If Needed, Train Custom MLP:
- Generate diverse DFT training data (~1000-10000 configs)
- Use modern architecture (MACE, NequIP)
- Apply active learning to fill gaps
- Validate extensively before production
Production:
- Monitor for instabilities during MD
- Compare key observables to experiment
- Report uncertainties in publications

10. Summary and Key Takeaways¶

What You Should Remember:¶

MLPs bridge the accuracy-speed gap
- DFT accuracy (~1 meV/atom)
- Near-classical speed (10-100× slower, not 1000×)
Descriptors encode local environments
- Must be invariant to translations, rotations, permutations
- Hand-crafted (symmetry functions) or learned (message passing)
Equivariant architectures are the state-of-the-art
- NequIP, MACE: 100× more data-efficient
- Foundation models enable zero-shot predictions
Training requires quality data and careful validation
- Diverse, relevant DFT calculations
- Active learning to fill gaps
- Validate with physical properties, not just MAE
Extrapolation is the Achilles’ heel
- MLPs can fail silently outside training distribution
- Always validate, never blindly trust
Foundation models are changing the field
- MACE-MP-0, CHGNet: pre-trained, general-purpose
- Good for screening, fine-tune for production

The Future:¶

Foundation models becoming standard starting points
Seamless DFT → MLP workflows
Uncertainty quantification becoming routine
MLPs + enhanced sampling for rare events
Multi-fidelity approaches (cheap + expensive data)

11. Resources and Further Reading¶

Software:¶

MLP Packages:

MACE - State-of-the-art equivariant GNN
NequIP - Equivariant neural networks
DeePMD-kit - Deep Potential, GPU-accelerated
n2p2 - Behler-Parrinello NNP
QUIP/GAP - Gaussian Approximation Potential

Foundation Models:

MACE-MP-0 - Materials Project model
CHGNet - Charge-informed GNN
M3GNet - Materials GNN

Training Data:

Materials Project - Inorganic materials
QM7/QM9 - Small organic molecules
ANI Dataset - Organic molecules

Key Papers:¶

Behler & Parrinello (2007) - “Generalized Neural-Network Representation of High-Dimensional Potential-Energy Surfaces” - The original NNP paper
Bartók et al. (2010) - “Gaussian Approximation Potentials” - Kernel methods with SOAP
Schütt et al. (2017) - “SchNet” - First message-passing NN for molecules
Batzner et al. (2022) - “E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials” - NequIP paper
Batatia et al. (2022) - “MACE: Higher Order Equivariant Message Passing Neural Networks” - Current SOTA

Reviews:¶

Unke et al. (2021) - “Machine Learning Force Fields” - Comprehensive review
Noé et al. (2020) - “Machine Learning for Molecular Simulation” - Broad overview
Friederich et al. (2021) - “Machine-learned potentials for next-generation matter simulations” - Perspective on future directions

12. Looking Ahead¶

In the tutorials, you will:

Use MACE-MP-0 to simulate a system without training
Compare MLP predictions to DFT reference calculations
Fine-tune a foundation model on custom data
Run production MD with an MLP potential in LAMMPS
Analyze the quality and limitations of MLP predictions

The era of machine learning potentials is here. While they require more expertise than classical force fields, they open doors to simulations that were previously impossible—reactive chemistry at large scale, accurate surface science, and materials discovery at quantum accuracy with classical speed.