Machine Learning Features/Descriptors, and Algorithms in Materials Science¶
1. Why Raw Structures Don't Work as Inputs¶
A crystal or molecule is just a list of atomic species and 3D coordinates. Algorithms need a fixed-format numerical object (a vector or a graph), and that object must obey three physical invariances:
- Translation invariance — shifting the whole structure in space shouldn't change the prediction.
- Rotation invariance — rotating the structure shouldn't change the prediction.
- Permutation invariance — relabeling atom 1 as atom 2 (same physical structure, different bookkeeping) shouldn't change the prediction.
A "descriptor" or "featurizer" is simply a function
that is built so these three symmetries are respected automatically, rather than something a model has to learn from data.
2. Structural Descriptors (Fingerprint Vectors)¶
2.1 The Coulomb Matrix (molecules / clusters)¶
For a system of \(N\) atoms, build an \(N \times N\) matrix \(M\):
Where each symbol comes from, physically:
- \(Z_i\): nuclear charge (atomic number) of atom \(i\) — heavier atoms get a bigger diagonal value.
- \(\lVert \mathbf{R}_i - \mathbf{R}_j \rVert\): Euclidean distance, \(\sqrt{(x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2}\).
- The off-diagonal form mirrors the classical Coulomb energy between two point charges, \(E \propto \frac{q_1 q_2}{r}\).
- The exponent \(2.4\) on the diagonal is empirical — it was fit so the self-term roughly reproduces the total electronic energy of a free atom.
Tiny worked example (2 atoms, H–F, \(Z_H=1\), \(Z_F=9\), bond length \(\approx 0.92\ \text{Å}\)):
Why permutation is a problem here: swapping atom order swaps rows/columns of \(M\), giving a different matrix for the same molecule. In practice this is fixed by sorting rows/columns by norm, or by using it only through permutation-invariant summaries (eigenvalues of \(M\), or a sorted Coulomb matrix).
Simple picture: \(M\) is a table where the diagonal says "how heavy is this atom" and every off-diagonal cell says "how strongly do these two atoms electrostatically feel each other."
2.2 Atom-Centered Symmetry Functions (ACSF)¶
Periodic crystals have no fixed atom count, so instead of one matrix for the whole system, ACSF builds a local fingerprint around each atom \(i\) from its neighbors \(j\).
Radial term (2-body, "how many atoms are at what distance"):
- \(R_{ij}\): distance from center atom \(i\) to neighbor \(j\).
- \(\eta\): controls how narrow the Gaussian probe is (large \(\eta\) = very local shell).
- \(R_s\): controls where that shell is centered (e.g. probe the 2nd nearest-neighbor shell).
- \(f_c\): smooth cutoff so atoms outside a physical radius \(R_c\) contribute nothing:
This cutoff smoothly decays to zero (and has zero slope at \(R_c\)), so as an atom drifts in or out of the cutoff sphere during a simulation, the descriptor — and any energy/force built from it — doesn't jump discontinuously.
Angular term (3-body, "at what angles do neighbor pairs sit"): this is what actually distinguishes, e.g., a tetrahedral site from an octahedral one at the same set of bond lengths.
- \(\theta_{ijk}\): the bond angle at atom \(i\) formed by neighbors \(j\) and \(k\).
- \(\lambda = \pm 1\): shifts which angles (\(0^\circ\) or \(180^\circ\)) are emphasized.
- \(\zeta\): sharpens the angular peak — larger \(\zeta\) means the function responds only to a narrow range of angles.
A full ACSF fingerprint for atom \(i\) is a vector built from many \(G^{\text{rad}}\) and \(G^{\text{ang}}\) terms, each with different \((\eta, R_s)\) or \((\eta,\zeta,\lambda)\) — effectively a bar-code of "how much stuff is nearby, at which distances, and at which angles."
3. Graph Representations (Crystal Graphs)¶
Represent the crystal as \(\mathcal{G} = (\mathcal{V}, \mathcal{E})\):
- Nodes \(\mathcal{V}\): atoms, each carrying a feature vector \(\mathbf{v}_i\) (atomic number, electronegativity, group/period, valence electron count, ...).
- Edges \(\mathcal{E}\): bonds, each carrying a feature vector \(\mathbf{e}_{ij}\) (bond length, sometimes bond angle, or a Gaussian-expanded version of \(R_{ij}\)).
Message passing update (one layer, schematically):
- \(\mathcal{N}(i)\): the set of bonded/nearby neighbors of atom \(i\).
- \(W_1, W_2, W_3\): learnable weight matrices.
- \(\sigma\): a nonlinearity (e.g. softplus, ReLU).
Each "message" \(\mathbf{m}_{ij}\) packages up what neighbor \(j\) tells atom \(i\); summing over neighbors and adding it back onto \(\mathbf{v}_i\) updates atom \(i\)'s representation using its local chemical environment. Stacking \(L\) such layers lets information from \(L\) bonds away reach a given atom — the graph analog of extending ACSF's cutoff radius. After several layers, node features are pooled (e.g. summed or averaged) into a single whole-crystal vector, which is then regressed against a target property.
Simple picture: the crystal is a subway map — stations (nodes) carry chemical identity, tracks (edges) carry bond distance, and "riding the tracks" for a few stops is how each atom learns about its neighborhood.
4. Connecting Descriptors to Supervised Learning¶
Once a structure is turned into \(\mathbf{x}\) (fingerprint vector or pooled graph embedding):
- Regression (continuous target, e.g. formation energy \(E_f\)): learn \(f(\mathbf{x}) \approx y\) using Random Forests, kernel ridge regression, neural networks, or — see below — Gaussian Process Regression.
- Classification (categorical target, e.g. space group): feed \(\mathbf{x}\) into an SVM or a softmax layer to get class probabilities \(p(y=c\mid\mathbf{x})\).
5. Gaussian Process Regression (GPR)¶
5.1 What "Gaussian" Actually Means¶
The ordinary 1D Gaussian (normal) distribution describes how likely a single number \(x\) is, given a mean \(\mu\) (center) and variance \(\sigma^2\) (spread):
It's the familiar bell curve: most probability mass sits near \(\mu\), and \(\sigma\) sets how fast it falls off.
This generalizes to a vector \(\mathbf{x} = (x_1,\dots,x_k)\) — a multivariate Gaussian:
Here \(\boldsymbol{\mu}\) is a vector of means (one per coordinate) and \(\Sigma\) is a \(k\times k\) covariance matrix: \(\Sigma_{ab} = \text{Cov}(x_a, x_b)\) tells you how strongly coordinate \(a\) and coordinate \(b\) move together. This one object — a covariance matrix — is the entire idea GPR is built on.
5.2 From Gaussian Vectors to Gaussian Processes¶
A Gaussian Process (GP) is the infinite-dimensional extension of the idea above: instead of a distribution over a finite vector, it's a distribution over functions. Any finite set of function values \(f(x_1), f(x_2), \dots, f(x_n)\) is assumed to be jointly Gaussian:
- \(m(\mathbf{x}) = \mathbb{E}[f(\mathbf{x})]\): the mean function (often just set to 0 after centering the data).
- \(k(\mathbf{x},\mathbf{x}') = \mathbb{E}\big[(f(\mathbf{x})-m(\mathbf{x}))(f(\mathbf{x}')-m(\mathbf{x}'))\big]\): the covariance/kernel function — this plays the exact same role as \(\Sigma_{ab}\) above, except now it's defined for any pair of inputs, not just a fixed finite list.
Intuition: a GP is a probability distribution over an infinite bag of candidate functions that could explain your data. Before seeing any data, all "reasonable" functions (as defined by the kernel) are equally in play. Every data point you observe rules out functions that don't pass near it, narrowing the bag.
5.3 The Kernel: Encoding "Nearby Inputs Should Give Similar Outputs"¶
The most common choice is the squared-exponential (RBF) kernel:
- \(\sigma_f^2\): signal variance — overall vertical scale of the function's variation.
- \(\ell\): length scale — how far apart two inputs can be before the model stops treating them as correlated. Small \(\ell\) → wiggly, locally-fit function. Large \(\ell\) → smooth, almost-linear function.
These are hyperparameters, usually tuned by maximizing the marginal likelihood of the observed data (not by cross-validating a loss the way you'd tune a neural net).
5.4 The Actual Prediction (Regression) Equations¶
Given training inputs \(X=\{\mathbf{x}_1,\dots,\mathbf{x}_n\}\) with (noisy) targets \(\mathbf{y}\), and a new test point \(\mathbf{x}_*\), GPR assumes:
where \(K_{ab}=k(\mathbf{x}_a,\mathbf{x}_b)\) is the \(n\times n\) training-training kernel matrix, \(\mathbf{k}_*\) is the vector \(k(\mathbf{x}_*,\mathbf{x}_a)\), \(k_{**}=k(\mathbf{x}_*,\mathbf{x}_*)\), and \(\sigma_n^2\) is assumed measurement noise. Conditioning this joint Gaussian on the observed \(\mathbf{y}\) (standard Gaussian conditioning identities) gives closed-form predictions:
Why this is elegant: the prediction is literally a weighted average of the observed \(\mathbf{y}\)'s, where the weights come from how "kernel-similar" \(\mathbf{x}_*\) is to each training point. And you get the variance — a rigorous, built-in error bar — for free, from the same formula. Far from any training point, \(\mathbf{k}_* \to 0\), so \(\bar f_* \to 0\) (falls back to the prior) and \(\mathbb{V}[f_*] \to k_{**}\) (uncertainty grows to its maximum).
5.5 Why GPR Matters for Materials/Catalysis Workflows¶
- Machine-learned interatomic potentials (GAP): Gaussian Approximation Potentials fit the local energy of an atomic environment (often described with an ACSF-like or SOAP descriptor as the kernel input) directly as a GP, giving forces/energies at a fraction of DFT cost, with an uncertainty estimate on every prediction.
- Active learning for expensive DFT/ab initio data: the posterior variance tells you exactly where your surrogate model is unreliable, so you can select the next structure to run in VASP/QE as the one with the highest predicted variance — squeezing maximum information out of every expensive calculation.
- Bayesian optimization for catalyst/composition screening: GPR is the surrogate model inside Bayesian optimization loops that search catalyst composition or adsorption-site space for optimal binding energy (e.g. sitting near the peak of a Sabatier volcano) using far fewer DFT evaluations than a grid search.
- Uncertainty-aware property prediction: unlike a plain neural network point-estimate, a GP tells you "trust this formation-energy prediction" vs. "this structure is far from anything in the training set, don't trust it" — valuable when screening large hypothetical structure databases.
One-line summary: a Gaussian Process treats "what function fits my data" itself as a random variable with a Gaussian distribution, and GPR is just the (closed-form) act of updating that distribution once you condition on observed points — giving you both a best-guess prediction and a mathematically honest confidence interval around it.
6. Quick Comparison¶
| Method | Input handled | Captures 3-body/angular info? | Native uncertainty? | Typical use |
|---|---|---|---|---|
| Coulomb Matrix | Molecules/clusters (fixed N) | No | No | Small-molecule property prediction |
| ACSF | Periodic/local environments | Yes (angular term) | No | Interatomic potentials, local energies |
| Crystal Graph (MPNN/CGCNN) | Arbitrary periodic structures | Yes (via edges + message passing) | No (unless made Bayesian) | Bulk property prediction (formation energy, band gap) |
| GPR (on top of any descriptor above) | Any fixed-length fingerprint | Depends on descriptor used | Yes, built-in | Surrogate potentials (GAP), active learning, Bayesian optimization |