ACE — Atomic Cluster Expansion descriptors#
Atomic Cluster Expansion (ACE) descriptors give every atom a fingerprint of its local environment. They are rotationally, translationally and permutationally invariant — exactly the properties an ML potential needs.
ACE was introduced by Drautz (2019). Several other popular descriptors are special cases of ACE:
descriptor |
corresponds to |
|---|---|
SOAP power spectrum |
ACE with \(\nu = 2\) |
SNAP / bispectrum |
ACE with \(\nu = 3\) |
Recipe#
Given an atom \(i\) and its neighbours \(j\) at displacement \(\vec{r}_{ij}\):
Project each bond onto a basis of radial \(R_n(r)\) and angular \(Y_l^m(\hat r)\) functions and sum to get the single-particle basis
$\(A_{i,nlm} = \sum_{j\in\mathrm{nbrs}(i)} R_n(r_{ij})\, Y_l^m(\hat r_{ij}).\)$Build rotationally invariant B-features by coupling products of \(A\)’s:
\(\nu = 1\): radial density \(B^{(1)}_{i,n} = A_{i,n00}\).
\(\nu = 2\): pair correlations (SOAP power spectrum) \(B^{(2)}_{i,n_1n_2l} = \sum_m A^*_{i,n_1lm}\,A_{i,n_2lm}\).
\(\nu = 3\): triplet correlations (bispectrum-like) using Clebsch–Gordan coupling.
Below we walk through the API and show that ACE features cleanly separate FCC, BCC and HCP environments — but only if you visualise them correctly. There is a small data-science gotcha that often hides the structural signal.
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from ase.build import bulk
import pyscal3 as pyscal
1. A first ACE calculation#
Build an FCC copper supercell, find neighbours and compute ACE descriptors with \(n_\max = 4\), \(l_\max = 6\), \(\nu_\max = 2\).
atoms = bulk("Cu", "fcc", cubic=True).repeat(3)
pyscal.find_neighbors(atoms, method="cutoff", cutoff=5.0)
result = pyscal.ace(atoms, nmax=4, lmax=6, nu_max=2)
print(f"{'block':<6} shape")
for key, val in result.items():
print(f"{key:<6} {val.shape}")
block shape
nu1 (108, 4)
nu2 (108, 70)
full (108, 74)
Each row of full is a fingerprint for one atom. The descriptor is split into
blocks by correlation order (\(\nu\)). Higher \(\nu\) captures higher-order angular
correlations but adds more features.
2. How the parameters control descriptor size#
ACE has three knobs:
nmax— radial resolution (number of Chebyshev radial basis functions).lmax— angular resolution (highest spherical-harmonic order).nu_max— correlation order (1 = density, 2 = pairs/SOAP, 3 = triplets).
Let’s sweep each one and watch the descriptor count grow.
small = bulk("Cu", "fcc", cubic=True).repeat(2)
pyscal.find_neighbors(small, method="cutoff", cutoff=5.0)
rows = []
for nmax in [2, 4, 6, 8]:
r = pyscal.ace(small, nmax=nmax, lmax=4, nu_max=2)
rows.append(("nmax", nmax, r['nu1'].shape[1], r['nu2'].shape[1], r['full'].shape[1]))
for lmax in [0, 2, 4, 6, 8]:
r = pyscal.ace(small, nmax=4, lmax=lmax, nu_max=2)
rows.append(("lmax", lmax, r['nu1'].shape[1], r['nu2'].shape[1], r['full'].shape[1]))
for nu in [1, 2, 3]:
r = pyscal.ace(small, nmax=4, lmax=4, nu_max=nu)
n2 = r['nu2'].shape[1] if 'nu2' in r else 0
n3 = r['nu3'].shape[1] if 'nu3' in r else 0
rows.append(("nu_max", nu, r['nu1'].shape[1], n2, n3))
print(f"{'sweep':<8} {'value':>5} {'nu1':>6} {'nu2':>6} {'nu3 / total':>14}")
for r in rows:
print(f"{r[0]:<8} {r[1]:>5} {r[2]:>6} {r[3]:>6} {r[4]:>14}")
sweep value nu1 nu2 nu3 / total
nmax 2 2 15 17
nmax 4 4 50 54
nmax 6 6 105 111
nmax 8 8 180 188
lmax 0 4 10 14
lmax 2 4 30 34
lmax 4 4 50 54
lmax 6 4 70 74
lmax 8 4 90 94
nu_max 1 4 0 0
nu_max 2 4 50 0
nu_max 3 4 50 110
Observations:
nu1only depends onnmax(it is just \(A_{i,n00}\)).nu2grows like \(n_\max^2 \,(l_\max+1)\) — angular detail and radial detail multiply.nu3grows fast; choose it only if you really need triplet information.
3. Inside the radial basis#
The radial functions are Chebyshev polynomials on \([r_\min, r_\mathrm{cut}]\)
multiplied by a smooth cutoff. Visualising them helps build intuition for what
nmax actually controls.
from pyscal3.descriptors import _ace_radial_basis
cutoff = 5.0
rs = np.linspace(0.5, cutoff, 400)
fig = go.Figure()
for n in range(6):
fig.add_trace(go.Scatter(x=rs, y=[_ace_radial_basis(n, r, cutoff) for r in rs],
mode="lines", name=f"R_{n}(r)"))
fig.update_layout(
title="ACE radial basis (Chebyshev × smooth cutoff)",
xaxis_title="r (\u00c5)", yaxis_title="R_n(r)",
height=380, template="simple_white",
)
fig
Each \(R_n\) has \(n\) zero-crossings. With \(n_\max = 6\) the radial direction can resolve about six independent shells before the cutoff.
4. Distinguishing FCC, BCC and HCP — the right way#
Now the critical part. Let’s compute ACE for three crystal structures and try to tell them apart from their fingerprints.
structures = {
"FCC (Cu)": bulk("Cu", "fcc", cubic=True).repeat(3),
"BCC (Fe)": bulk("Fe", "bcc", cubic=True).repeat(3),
"HCP (Mg)": bulk("Mg", "hcp").repeat((3, 3, 2)),
}
per_atom = {}
for name, atoms in structures.items():
pyscal.find_neighbors(atoms, method="cutoff", cutoff=5.0)
res = pyscal.ace(atoms, nmax=4, lmax=6, nu_max=2, normalize=False)
per_atom[name] = res['full']
print(f"{name}: {res['full'].shape[0]} atoms × {res['full'].shape[1]} features")
FCC (Cu): 108 atoms × 74 features
BCC (Fe): 54 atoms × 74 features
HCP (Mg): 36 atoms × 74 features
The naive plot — and why it misleads#
If you simply average the descriptor across atoms and bar-plot the raw values, FCC and BCC look almost indistinguishable. The reason is purely visual: a couple of large radial components (the average neighbour count) dominate the y-axis and dwarf the angular components that actually differ between structures.
fig = make_subplots(rows=1, cols=3, subplot_titles=list(per_atom.keys()),
shared_yaxes=True)
for col, (name, X) in enumerate(per_atom.items(), start=1):
fp = X.mean(axis=0)
fig.add_trace(go.Bar(x=np.arange(len(fp)), y=fp, name=name,
marker_line_width=0, showlegend=False), row=1, col=col)
fig.update_layout(height=320, title="Raw mean ACE fingerprint (misleading)",
template="simple_white")
fig.update_xaxes(title="feature index")
fig.update_yaxes(title="value", col=1)
fig
Compare just the correlations:
names = list(per_atom)
raw = np.stack([X.mean(axis=0) for X in per_atom.values()])
C = np.corrcoef(raw)
print("Pearson correlation between raw mean fingerprints:")
print(f"{'':10s}" + "".join(f"{n:>12s}" for n in names))
for i, n in enumerate(names):
print(f"{n:10s}" + "".join(f"{C[i,j]:>12.4f}" for j in range(3)))
Pearson correlation between raw mean fingerprints:
FCC (Cu) BCC (Fe) HCP (Mg)
FCC (Cu) 1.0000 0.9994 0.6762
BCC (Fe) 0.9994 1.0000 0.6717
HCP (Mg) 0.6762 0.6717 1.0000
FCC vs BCC reads as 0.999 — visually they look identical. That is an artefact of the magnitude scale, not of the descriptor. Each feature has its own characteristic scale, so we have to put them on equal footing first.
The fix — z-score across all atoms#
Stack the per-atom fingerprints from every structure together, subtract the global mean and divide by the global standard deviation. Now every feature contributes equally, and structural differences pop out.
labels = np.concatenate([np.full(len(X), name) for name, X in per_atom.items()])
X_all = np.vstack(list(per_atom.values()))
X_z = (X_all - X_all.mean(axis=0)) / (X_all.std(axis=0) + 1e-12)
fig = make_subplots(rows=1, cols=3, subplot_titles=names, shared_yaxes=True)
for col, name in enumerate(names, start=1):
fp = X_z[labels == name].mean(axis=0)
fig.add_trace(go.Bar(x=np.arange(len(fp)), y=fp, marker_line_width=0,
showlegend=False), row=1, col=col)
fig.update_layout(height=320,
title="Standardised mean ACE fingerprint — clear differences",
template="simple_white")
fig.update_xaxes(title="feature index")
fig.update_yaxes(title="z-score", col=1)
fig
Now FCC, BCC and HCP have visibly different patterns — and importantly, the angular features (high indices) light up too.
A 2-D map of every atom#
PCA on the standardised features projects every atom onto two coordinates. All atoms inside a perfect crystal collapse to a single point, and the three structures land at well-separated locations.
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
Z = pca.fit_transform(X_z)
fig = go.Figure()
for name in names:
pts = Z[labels == name]
fig.add_trace(go.Scatter(x=pts[:, 0], y=pts[:, 1], mode="markers",
name=name, marker=dict(size=10, opacity=0.7)))
fig.update_layout(
title=f"PCA of standardised ACE fingerprints "
f"(explained variance: {pca.explained_variance_ratio_.sum():.0%})",
xaxis_title="PC1", yaxis_title="PC2",
height=420, template="simple_white",
)
fig
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[9], line 1
----> 1 from sklearn.decomposition import PCA
2
3 pca = PCA(n_components=2)
4 Z = pca.fit_transform(X_z)
ModuleNotFoundError: No module named 'sklearn'
Each cluster is a single dot because every atom in a perfect crystal sees the exact same environment. PCA captures essentially 100 % of the variance with two components.
5. Distinguishing crystals from a melt#
A more practical use-case: detecting which atoms in a snapshot are crystalline and which are liquid-like. We perturb the FCC and BCC supercells with random displacements (a cheap proxy for thermal noise) and pretend a third class is a fully randomised “liquid”.
rng = np.random.default_rng(0)
def perturb(atoms, sigma):
a = atoms.copy()
a.positions += rng.normal(scale=sigma, size=a.positions.shape)
return a
fcc_T = perturb(bulk("Cu", "fcc", cubic=True).repeat(3), sigma=0.10)
bcc_T = perturb(bulk("Fe", "bcc", cubic=True).repeat(3), sigma=0.10)
liq = perturb(bulk("Cu", "fcc", cubic=True).repeat(3), sigma=0.6)
snapshots = {"FCC + noise": fcc_T, "BCC + noise": bcc_T, "liquid-like": liq}
X_list, y_list = [], []
for name, atoms in snapshots.items():
pyscal.find_neighbors(atoms, method="cutoff", cutoff=5.0)
res = pyscal.ace(atoms, nmax=4, lmax=6, nu_max=2, normalize=False)
X_list.append(res['full'])
y_list.append(np.full(len(atoms), name))
X = np.vstack(X_list)
y = np.concatenate(y_list)
X = (X - X.mean(axis=0)) / (X.std(axis=0) + 1e-12)
Z = PCA(n_components=2).fit_transform(X)
fig = go.Figure()
for name in snapshots:
pts = Z[y == name]
fig.add_trace(go.Scatter(x=pts[:, 0], y=pts[:, 1], mode="markers",
name=name, marker=dict(size=8, opacity=0.6)))
fig.update_layout(title="PCA of ACE fingerprints — noisy crystals vs liquid",
xaxis_title="PC1", yaxis_title="PC2",
height=420, template="simple_white")
fig
Even with 0.1 Å of noise the FCC and BCC clouds remain compact and well separated, while the liquid-like atoms scatter across a much wider region. This is exactly the kind of signal an ML potential or a structure classifier exploits.
6. Triplet correlations (\(\nu = 3\))#
Setting nu_max=3 adds bispectrum-like features. They are more expensive but
carry information that pair correlations cannot — particularly the handedness
and triangular geometry of three-body arrangements.
atoms = bulk("Cu", "fcc", cubic=True).repeat(3)
pyscal.find_neighbors(atoms, method="cutoff", cutoff=5.0)
res = pyscal.ace(atoms, nmax=4, lmax=4, nu_max=3, normalize=False)
fig = make_subplots(rows=1, cols=3, subplot_titles=["\u03bd = 1", "\u03bd = 2", "\u03bd = 3"])
for col, key in enumerate(["nu1", "nu2", "nu3"], start=1):
data = res[key].ravel()
fig.add_trace(go.Histogram(x=data, nbinsx=40, showlegend=False), row=1, col=col)
fig.update_layout(height=320, template="simple_white",
title="Distribution of ACE feature values per correlation order")
fig.update_xaxes(title="value")
fig.update_yaxes(title="count", col=1)
fig
7. Where results are stored#
After calling pyscal.ace, the descriptors and the parameters used are
attached to the Atoms object so they travel with the data.
print("atoms.arrays['pyscal_ace'] :", atoms.arrays['pyscal_ace'].shape)
print("atoms.info['pyscal_ace_params'] :", atoms.info['pyscal_ace_params'])
atoms.arrays['pyscal_ace'] : (108, 164)
atoms.info['pyscal_ace_params'] : {'nmax': 4, 'lmax': 4, 'nu_max': 3, 'cutoff': 5.0}
Take-aways#
ACE gives a complete, invariant fingerprint of an atom’s environment.
nmax,lmaxandnu_maxtrade resolution against descriptor size.Always standardise ACE features (or any high-dimensional descriptor) before plotting or computing similarities — raw means are dominated by a few large radial components and hide structural differences.
Once standardised, simple PCA cleanly separates FCC, BCC, HCP and liquid environments.
References#
Drautz, R. (2019). Atomic cluster expansion for accurate and transferable interatomic potentials. Phys. Rev. B 99, 014104.
Dusson et al. (2022). Atomic cluster expansion: completeness, efficiency and stability. J. Comput. Phys.
Bartók et al. (2010). Gaussian approximation potentials. Phys. Rev. Lett. 104, 136403.