Steinhardt Bond-Order Parameters

Steinhardt parameters \(q_l\) quantify the local symmetry around each atom using spherical harmonics. They are widely used to distinguish crystal structures:

Structure	\(q_4\)	\(q_6\)
FCC	0.19	0.57
BCC	0.04	0.51
HCP	0.10	0.48

import pyscal
from ase.build import bulk
import numpy as np

Basic Usage

# Create FCC Cu
atoms = bulk("Cu", "fcc", cubic=True).repeat(4)

# Find neighbors first
pyscal.find_neighbors(atoms, method="cutoff", cutoff=0)

# Compute q4 and q6
q = pyscal.steinhardt_parameter(atoms, l=[4, 6])

print(f"q4: {q[0].mean():.4f} +/- {q[0].std():.6f}")
print(f"q6: {q[1].mean():.4f} +/- {q[1].std():.6f}")

q4: 0.1909 +/- 0.000000
q6: 0.5745 +/- 0.000000

Comparing Crystal Structures

from pyscal.structures import make_crystal

structures = {
    "FCC": make_crystal("fcc", lattice_constant=3.6, repetitions=(4, 4, 4)),
    "BCC": make_crystal("bcc", lattice_constant=2.87, repetitions=(4, 4, 4)),
    "HCP": make_crystal("hcp", lattice_constant=3.2, repetitions=(4, 4, 4), ca_ratio=1.633),
}

print(f"{'Structure':<10} {'q4':>8} {'q6':>8}")
print("-" * 28)
for name, s in structures.items():
    pyscal.find_neighbors(s, method="cutoff", cutoff=0)
    q = pyscal.steinhardt_parameter(s, l=[4, 6])
    print(f"{name:<10} {q[0].mean():>8.4f} {q[1].mean():>8.4f}")

Structure        q4       q6
----------------------------
FCC          0.1909   0.5745
BCC          0.0364   0.5107
HCP          0.0972   0.4848

Averaged Steinhardt Parameters

Averaged parameters \(\bar{q}_l\) smooth out thermal noise by averaging the complex spherical harmonic coefficients over the neighbor shell. Set averaged=True.

atoms = bulk("Cu", "fcc", cubic=True).repeat(4)
pyscal.find_neighbors(atoms, method="cutoff", cutoff=0)

q_avg = pyscal.steinhardt_parameter(atoms, l=[4, 6], averaged=True)

print(f"Averaged q4: {q_avg[0].mean():.4f}")
print(f"Averaged q6: {q_avg[1].mean():.4f}")

Averaged q4: 0.1909
Averaged q6: 0.5745

Higher-Order Parameters

You can compute any \(l\) value. Higher values capture finer angular details.

atoms = bulk("Cu", "fcc", cubic=True).repeat(4)
pyscal.find_neighbors(atoms, method="cutoff", cutoff=0)

q = pyscal.steinhardt_parameter(atoms, l=[2, 4, 6, 8, 10, 12])

for i, l_val in enumerate([2, 4, 6, 8, 10, 12]):
    print(f"q{l_val:<3}: {q[i].mean():.4f}")

q2  : 0.0000
q4  : 0.1909
q6  : 0.5745
q8  : 0.4039
q10 : 0.0129
q12 : 0.6001

Accessing Stored Results

Results are stored in atoms.arrays with keys pyscal_q<l> (or pyscal_q<l>_avg for averaged).

# Direct access to per-atom arrays
q6_values = atoms.arrays["pyscal_q6"]
print(f"q6 per atom (first 5): {q6_values[:5]}")
print(f"Shape: {q6_values.shape}")

q6 per atom (first 5): [0.57452426 0.57452426 0.57452426 0.57452426 0.57452426]
Shape: (256,)