Minkowski structure metrics#
Standard Steinhardt \(q_l\) counts every neighbour equally — but “who is a neighbour” depends on a cutoff radius you have to pick. Two atoms sitting right at the cutoff jump in or out of the count, which makes \(q_l\) jittery for noisy or thermally agitated configurations.
Minkowski structure metrics (Mickel et al., 2013) fix this by weighting each neighbour by the area of the Voronoi face they share. Distant neighbours have small faces and contribute little; close neighbours dominate. Crucially, no cutoff is needed and the result varies smoothly under thermal motion.
where \(A_{ij}\) is the Voronoi face shared between \(i\) and \(j\).
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. Perfect crystals#
On an ideal lattice every atom has the same Voronoi cell, so every atom gets the same Minkowski \(q_l\). The values are characteristic of the structure.
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():
q4, q6 = pyscal.minkowski_parameter(atoms, l=[4, 6])
rows.append((name, q4[0], q6[0]))
print(f"{'':5s} {'q4_Mink':>10} {'q6_Mink':>10}")
for name, q4, q6 in rows:
print(f'{name:5s} {q4:>10.4f} {q6:>10.4f}')
q4_Mink q6_Mink
FCC 0.1909 0.5745
BCC 0.2240 0.5669
HCP 0.0972 0.4848
2. Standard \(q_l\) vs Minkowski \(q_l\) — the cutoff problem#
Compute the standard \(q_6\) for FCC across a range of cutoff radii. As the cutoff sweeps from “first shell only” to “first + second shells”, the value jumps. The Minkowski version stays put.
fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4))
# Minkowski q6 (no cutoff to choose)
[q6_mink] = pyscal.minkowski_parameter(fcc, l=6)
q6_mink_value = float(q6_mink[0])
cutoffs = np.linspace(2.6, 4.5, 30)
q6_std = []
for c in cutoffs:
a = fcc.copy()
pyscal.find_neighbors(a, method='cutoff', cutoff=c)
[q6] = pyscal.steinhardt_parameter(a, l=6)
q6_std.append(q6.mean())
fig = go.Figure()
fig.add_trace(go.Scatter(x=cutoffs, y=q6_std, mode='lines+markers',
name='standard q6 (cutoff dependent)'))
fig.add_hline(y=q6_mink_value, line_dash='dash', line_color='crimson',
annotation_text=f'Minkowski q6 = {q6_mink_value:.3f}',
annotation_position='bottom right')
fig.update_layout(xaxis_title='cutoff (\u00c5)', yaxis_title='q6',
title='Standard q6 jumps with cutoff; Minkowski q6 is constant',
height=380, template='simple_white')
fig
Notice the abrupt jump near 3.4 Å — that’s the second neighbour shell entering the count. The Minkowski curve is just a flat line because face areas decay smoothly as new neighbours appear.
3. Behaviour under thermal noise#
Add increasing Gaussian displacements to FCC and look at how the per-atom \(q_6^{\text{Mink}}\) distribution broadens. Even at strong noise the central value stays near the perfect-crystal answer.
noises = [0.0, 0.02, 0.05, 0.10, 0.20]
fig = go.Figure()
stats = []
for s in noises:
a = make_crystal('fcc', lattice_constant=4.05, repetitions=(5, 5, 5), noise=s)
[q6] = pyscal.minkowski_parameter(a, l=6)
fig.add_trace(go.Histogram(x=q6, nbinsx=40, name=f'noise={s:.2f}',
opacity=0.55, histnorm='probability density'))
stats.append((s, float(q6.mean()), float(q6.std())))
fig.update_layout(barmode='overlay', xaxis_title='q6_Mink',
yaxis_title='density', template='simple_white',
title='Distribution of per-atom q6_Mink with thermal noise',
height=400)
fig
print(f"{'noise':>7} {'mean q6':>10} {'std':>10}")
for s, m, sd in stats:
print(f'{s:>7.2f} {m:>10.5f} {sd:>10.5f}')
noise mean q6 std
0.00 0.57452 0.00000
0.02 0.57021 0.00267
0.05 0.54965 0.01552
0.10 0.47484 0.05237
0.20 0.39577 0.07067
4. The face-area exponent voroexp#
voroexp = α raises each face area to the power \(\alpha\) before
normalisation:
\(\alpha = 0\): every Voronoi neighbour gets equal weight (recovers a cutoff-free Steinhardt \(q_l\)).
\(\alpha = 1\): standard Minkowski metric (default).
\(\alpha > 1\): emphasises the largest faces even more.
fcc = make_crystal('fcc', lattice_constant=4.05, repetitions=(4, 4, 4))
alphas = [0, 1, 2, 3, 4]
vals = [float(pyscal.minkowski_parameter(fcc, l=6, voroexp=a)[0][0]) for a in alphas]
fig = go.Figure(go.Scatter(x=alphas, y=vals, mode='lines+markers',
marker=dict(size=10)))
fig.update_layout(xaxis_title='voroexp \u03b1',
yaxis_title='q6_Mink',
title='Sensitivity to face-area exponent',
height=350, template='simple_white')
fig
5. Putting FCC, BCC, HCP on a single map#
The \((q_4^{\text{Mink}}, q_6^{\text{Mink}})\) plane is a clean fingerprint for crystal classification — and because it is cutoff-free, it transfers directly to noisy MD snapshots.
fig = go.Figure()
for name, lc, reps in [('FCC', 4.05, (4,4,4)),
('BCC', 2.87, (4,4,4)),
('HCP', 3.21, (4,4,4))]:
for sigma in [0.0, 0.10]:
a = make_crystal(name.lower(), lattice_constant=lc,
repetitions=reps, noise=sigma)
q4, q6 = pyscal.minkowski_parameter(a, l=[4, 6])
fig.add_trace(go.Scatter(x=q4, y=q6, mode='markers',
name=f'{name} (\u03c3={sigma:.2f})',
marker=dict(size=6, opacity=0.5)))
fig.update_layout(xaxis_title='q4_Mink', yaxis_title='q6_Mink',
title='Crystal-structure map with thermal noise',
template='simple_white', height=460)
fig
Take-aways#
Minkowski \(q_l\) is a cutoff-free alternative to Steinhardt \(q_l\): Voronoi face areas provide an automatic, smooth weighting.
Distributions are tighter under thermal noise than the standard Steinhardt \(q_l\).
The trade-off is cost: a Voronoi tessellation per atom, slower than a cutoff-based neighbour search.
References#
W. Mickel, S. C. Kapfer, G. E. Schröder-Turk and K. Mecke, Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter, J. Chem. Phys. 138, 044501 (2013).
P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Bond-orientational order in liquids and glasses, Phys. Rev. B 28, 784 (1983).