Wigner \(W_l\) parameters#
If you have used the Steinhardt \(q_l\) parameters before, the Wigner \(W_l\) are their natural cousins: where \(q_l\) measures the magnitude of bond orientational order, \(W_l\) measures a third-order structural invariant that often distinguishes structures \(q_l\) cannot.
Formally,
where the bracketed object is a Wigner 3j symbol. The normalised version \(\hat W_l\) divides by \(\big(\sum_m |q_{lm}|^2\big)^{3/2}\), so it depends only on the shape of the local environment, not on its overall order parameter magnitude.
Why we care: \(\hat W_6\) has opposite signs for FCC (\(-0.013\)) and BCC (\(+0.013\)). That single sign flip is enough to tell the two apart at a glance — something \(q_6\) alone cannot do because both lattices give very similar \(q_6\) values.
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import pyscal3 as pyscal
from pyscal3.structures import make_crystal
1. A first calculation#
Build a perfect FCC crystal, find neighbours and ask for \(\hat W_4\) and \(\hat W_6\).
fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4))
pyscal.find_neighbors(fcc, method='cutoff', cutoff=0) # adaptive cutoff
what4, what6 = pyscal.wigner_w_parameter(fcc, l=[4, 6])
print(f'FCC perfect crystal:')
print(f' W-hat_4 = {what4.mean():+.5f} (reference: -0.15932)')
print(f' W-hat_6 = {what6.mean():+.5f} (reference: -0.01316)')
FCC perfect crystal:
W-hat_4 = -0.15932 (reference: -0.15932)
W-hat_6 = -0.01316 (reference: -0.01316)
Both numbers match the textbook reference values from Steinhardt, Nelson & Ronchetti (1983) to within rounding.
2. The fingerprint of three crystals#
Now compute the same numbers for FCC, BCC and HCP and compare against the ordinary \(q_6\). Watch what \(q_6\) vs \(\hat W_6\) does for FCC vs BCC.
structures = {
'FCC': make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4)),
'BCC': make_crystal('bcc', lattice_constant=2.87, repetitions=(4, 4, 4)),
'HCP': make_crystal('hcp', lattice_constant=3.21, repetitions=(4, 4, 4)),
}
rows = []
for name, atoms in structures.items():
pyscal.find_neighbors(atoms, method='cutoff', cutoff=0)
w4, w6 = pyscal.wigner_w_parameter(atoms, l=[4, 6])
[q6] = pyscal.steinhardt_parameter(atoms, l=6)
rows.append((name, w4.mean(), w6.mean(), q6.mean()))
print(f"{'':5s} {'W-hat_4':>10} {'W-hat_6':>10} {'q_6':>10}")
for name, w4, w6, q6 in rows:
print(f'{name:5s} {w4:>+10.5f} {w6:>+10.5f} {q6:>10.4f}')
W-hat_4 W-hat_6 q_6
FCC -0.15932 -0.01316 0.5745
BCC +0.15932 +0.01316 0.5107
HCP +0.13410 -0.01244 0.4848
Read it like this:
\(q_6\): FCC and BCC are close (0.575 vs 0.511), HCP a bit lower. Hard to tell FCC and BCC apart from \(q_6\) alone.
\(\hat W_4\): opposite signs for FCC and BCC.
\(\hat W_6\): opposite signs again — and HCP is also negative but with a very different \(\hat W_4\).
So the \((\hat W_4, \hat W_6)\) pair is a clean two-coordinate fingerprint.
fig = go.Figure()
for name, w4, w6, _ in rows:
fig.add_trace(go.Scatter(x=[w4], y=[w6], mode='markers+text',
text=[name], textposition='top center',
marker=dict(size=18), name=name))
fig.add_hline(y=0, line_dash='dot', line_color='grey')
fig.add_vline(x=0, line_dash='dot', line_color='grey')
fig.update_layout(title='(W-hat_4, W-hat_6) fingerprint of perfect crystals',
xaxis_title='W-hat_4', yaxis_title='W-hat_6',
height=420, template='simple_white', showlegend=False)
fig
3. Why odd \(l\) vanishes#
The Wigner 3j symbol \(\begin{pmatrix} l & l & l \\ m_1 & m_2 & m_3 \end{pmatrix}\) is non-zero only when \(3l\) is even, i.e. when \(l\) itself is even. So \(W_l\) identically vanishes for odd \(l\) — a useful sanity check.
fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(3, 3, 3))
pyscal.find_neighbors(fcc, method='cutoff', cutoff=0)
for l in range(3, 9):
[w] = pyscal.wigner_w_parameter(fcc, l=l)
flag = ' (zero — odd l)' if l % 2 else ''
print(f'l={l}: W-hat_{l} = {w.mean():+.5f}{flag}')
l=3: W-hat_3 = +0.00000 (zero — odd l)
l=4: W-hat_4 = -0.15932
l=5: W-hat_5 = +0.00000 (zero — odd l)
l=6: W-hat_6 = -0.01316
l=7: W-hat_7 = +0.00000 (zero — odd l)
l=8: W-hat_8 = +0.05845
4. Survival under thermal noise#
Real MD snapshots are noisy. Two tools help:
The averaged \(\bar W_l\) (Lechner & Dellago, 2008) replaces \(q_{lm}\) by the average over an atom and its neighbours — this dramatically sharpens distributions.
A 2-D scatter on \((\hat W_4, \hat W_6)\) shows the noise as cloud width.
Add 0.05 Å of Gaussian noise to perfect FCC and BCC supercells and look at every atom’s signature.
fcc_n = make_crystal('fcc', lattice_constant=4.05, repetitions=(5, 5, 5), noise=0.05)
bcc_n = make_crystal('bcc', lattice_constant=2.87, repetitions=(5, 5, 5), noise=0.05)
data = {}
for name, atoms in [('FCC', fcc_n), ('BCC', bcc_n)]:
pyscal.find_neighbors(atoms, method='cutoff', cutoff=0)
w4_per, w6_per = pyscal.wigner_w_parameter(atoms, l=[4, 6])
w4_avg, w6_avg = pyscal.wigner_w_parameter(atoms, l=[4, 6], averaged=True)
data[name] = (w4_per, w6_per, w4_avg, w6_avg)
fig = make_subplots(rows=1, cols=2,
subplot_titles=('per-atom W-hat',
'neighbour-averaged \u00afW-hat'),
shared_xaxes=False, shared_yaxes=False)
for name, (w4p, w6p, w4a, w6a) in data.items():
fig.add_trace(go.Scatter(x=w4p, y=w6p, mode='markers',
marker=dict(size=4, opacity=0.4), name=name,
legendgroup=name), row=1, col=1)
fig.add_trace(go.Scatter(x=w4a, y=w6a, mode='markers',
marker=dict(size=4, opacity=0.4), name=name,
legendgroup=name, showlegend=False), row=1, col=2)
fig.update_xaxes(title='W-hat_4')
fig.update_yaxes(title='W-hat_6')
fig.add_hline(y=0, line_dash='dot', line_color='grey', row=1, col=1)
fig.add_vline(x=0, line_dash='dot', line_color='grey', row=1, col=1)
fig.add_hline(y=0, line_dash='dot', line_color='grey', row=1, col=2)
fig.add_vline(x=0, line_dash='dot', line_color='grey', row=1, col=2)
fig.update_layout(height=420, template='simple_white',
title='Noisy FCC vs BCC in the (W-hat_4, W-hat_6) plane')
fig
Per-atom values smear into broad clouds, but the FCC and BCC clusters sit in opposite quadrants and remain separable. The averaged variant (\(\bar W\)) collapses each cloud onto a tight blob — much easier for any downstream classifier or threshold rule.
5. Normalised vs raw \(W_l\)#
By default wigner_w_parameter returns the dimensionless \(\hat W_l\).
Pass normalized=False to get the raw \(W_l\) (useful when you want a
physically sized quantity rather than an angle-only invariant).
fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(3, 3, 3))
pyscal.find_neighbors(fcc, method='cutoff', cutoff=0)
[w_hat] = pyscal.wigner_w_parameter(fcc, l=6, normalized=True)
[w_raw] = pyscal.wigner_w_parameter(fcc, l=6, normalized=False)
print(f'normalised W-hat_6 = {w_hat.mean():+.5f}')
print(f'raw W_6 = {w_raw.mean():+.5e}')
normalised W-hat_6 = -0.01316
raw W_6 = -2.62604e-03
6. Where it is stored#
Each call writes into atoms.arrays under a pyscal_ prefix so the data
follows the structure around.
fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(3, 3, 3))
pyscal.find_neighbors(fcc, method='cutoff', cutoff=0)
pyscal.wigner_w_parameter(fcc, l=[4, 6], averaged=True)
for key in sorted(fcc.arrays):
if 'w' in key.lower() and key.startswith('pyscal_'):
print(f' {key:25s} shape={fcc.arrays[key].shape} mean={fcc.arrays[key].mean():+.5f}')
pyscal_avg_w4 shape=(108,) mean=-0.00067
pyscal_avg_w6 shape=(108,) mean=-0.00263
pyscal_avg_what4 shape=(108,) mean=-0.15932
pyscal_avg_what6 shape=(108,) mean=-0.01316
pyscal_neighborweight shape=(108, 12) mean=+1.00000
Take-aways#
\(\hat W_l\) is a third-order invariant: same input data as \(q_l\), but carries the sign and shape information that \(q_l\) throws away.
The pair \((\hat W_4, \hat W_6)\) is a clean fingerprint that separates FCC, BCC, HCP and icosahedral environments.
For finite-temperature snapshots, use the averaged variant — same call, with
averaged=True.
References#
Steinhardt, Nelson & Ronchetti, Bond-orientational order in liquids and glasses, Phys. Rev. B 28, 784 (1983).
Lechner & Dellago, Accurate determination of crystal structures based on averaged local bond order parameters, J. Chem. Phys. 129, 114707 (2008).