freud.environment.LocalDescriptors: Steinhardt Order Parameters from Scratch

The freud.environment module analyzes the local environments of particles. The freud.environment.LocalDescriptors class is a useful tool for analyzing identifying crystal structures in a rotationally invariant manner using local particle environments. The primary purpose of this class is to compute spherical harmonics between neighboring particles in a way that orients particles correctly relative to their local environment, ensuring that global orientational shifts do not change the output.

import freud
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

Computing Spherical Harmonics

To demonstrate the basic application of the class, let’s compute the spherical harmonics between neighboring particles. For simplicity, we consider points on a simple cubic lattice.

uc =
box, points = uc.generate_system(5)
system = freud.AABBQuery(box, points)
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection="3d")
ax.set_title("Simple cubic crystal", fontsize=16)

Now, let’s use the class to compute an array of spherical harmonics for the system. The harmonics are computed for each bond, where a bond is defined by a pair of particles that are determined to lie within each others’ nearest neighbor shells based on a standard neighbor list search. The number of bonds and spherical harmonics to calculate is configurable.

num_neighbors = 6
l_max = 12

nlist = system.query(
    points, {"num_neighbors": num_neighbors, "exclude_ii": True}
ld = freud.environment.LocalDescriptors(l_max, mode="global")
ld.compute(system, neighbors=nlist);

Accessing the Data

The resulting spherical harmonic array has a shape corresponding to the number of neighbors. We can now extract the spherical harmonics corresponding to a particular \((l, m)\) pair using the ordering used by the LocalDescriptors class: increasing values of \(l\), and for each \(l\), the nonnegative \(m\) values followed by the negative values.

sph_raw = np.mean(ld.sph, axis=0)
count = 0
sph = np.zeros((l_max + 1, l_max + 1), dtype=np.complex128)
for l in range(l_max + 1):
    for m in range(l + 1):
        sph[l, m] = sph_raw[count]
        count += 1
    for m in range(-l, 0):
        sph[l, m] = sph_raw[count]
        count += 1

Using Spherical Harmonics to Compute Steinhardt Order Parameters

The raw per bond spherical harmonics are not typically useful quantities on their own. However, they can be used to perform sophisticated crystal structure analyses with different methods; for example, the pythia library uses machine learning to find patterns in the spherical harmonics computed by this class. In this notebook, we’ll use the quantities for a more classical application: the computation of Steinhardt order parameters. The order parameters \(q_l\) provide a rotationally invariant measure of the system that can for some structures, provide a unique identifying fingerprint. They are a particularly useful measure for various simple cubic structures such as structures with underlying simple cubic, BCC, or FCC lattices. The freud library actually provides additional classes to efficiently calculate these order parameters directly, but they also provide a reasonable demonstration here.

For more information on Steinhardt order parameters, see the original paper or the freud.order.Steinhardt documentation.

def get_ql(num_particles, descriptors, nlist, weighted=False):
    """Given a set of points and a LocalDescriptors object (and the
    underlying NeighborList), compute the per-particle Steinhardt ql
    order parameter for all :math:`l` values up to the maximum quantum
    number used in the computation of the descriptors."""
    qbar_lm = np.zeros((num_particles, descriptors.sph.shape[1]), dtype=np.complex128)
    for i in range(num_particles):
        indices = nlist.query_point_indices == i
        Ylms = descriptors.sph[indices, :]
        if weighted:
            weights = nlist.weights[indices, np.newaxis]
            weights /= np.sum(weights)
            num_neighbors = 1
            weights = np.ones_like(Ylms)
            num_neighbors = descriptors.sph.shape[0] / num_particles
        qbar_lm[i, :] = np.sum(Ylms * weights, axis=0) / num_neighbors

    ql = np.zeros((qbar_lm.shape[0], descriptors.l_max + 1))
    for i in range(ql.shape[0]):
        for l in range(ql.shape[1]):
            for k in range(l**2, (l + 1) ** 2):
                ql[i, l] += np.absolute(qbar_lm[i, k]) ** 2
            ql[i, l] = np.sqrt(4 * np.pi / (2 * l + 1) * ql[i, l])

    return ql

ld_ql = get_ql(len(points), ld, nlist)

Since freud provides the ability to calculate these parameter as well, we can directly check that our answers are correct. Note: More information on the ``Steinhardt`` class can be found in the documentation or in the ``Steinhardt`` example.

L = 6
steinhardt = freud.order.Steinhardt(l=L)
steinhardt.compute(system, neighbors=nlist)
if np.allclose(steinhardt.ql, ld_ql[:, L]):
    print("Our manual calculation matches the Steinhardt class!")
Our manual calculation matches the Steinhardt class!

For a brief demonstration of why the Steinhardt order parameters can be useful, let’s look at the result of thermalizing our points and recomputing this measure.

sigmas = [0.03, 0.05, 0.1]
systems = []
nlists = []
for sigma in sigmas:
    box, points = uc.generate_system(5, sigma_noise=sigma)
    system = freud.AABBQuery(box, points)
            points, {"num_neighbors": num_neighbors, "exclude_ii": True}
fig = plt.figure(figsize=(14, 6))
axes = []
for i, v in enumerate(sigmas):
    ax = fig.add_subplot(int(f"1{len(sigmas)}{i+1}"), projection="3d")
    ax.set_title(rf"$\sigma$ = {v}", fontsize=16)

If we recompute the Steinhardt OP for each of these data sets, we see that adding noise has the effect of smoothing the order parameter such that the peak we observed for the perfect crystal is no longer observable.

ld_qls = []
for i, sigma in enumerate(sigmas):
    ld = freud.environment.LocalDescriptors(l_max, mode="global")
    ld.compute(systems[i], neighbors=nlists[i])
    ld_qls.append(get_ql(len(systems[i].points), ld, nlists[i]))
fig, ax = plt.subplots()
for i, ld_ql in enumerate(ld_qls):
    lim_out = ax.hist(ld_ql[:, L], label=rf"$\sigma$ = {sigmas[i]}", density=True)
    if i == 0:
        # Can choose any element, all are identical in the reference case
        ax.vlines(ld_ql[:, L][0], 0, np.max(lim_out[0]), label="Reference")
ax.set_title(f"Histogram of $q_{L}$ values", fontsize=16)
ax.set_ylabel("Frequency", fontsize=14)
ax.set_xlabel(f"$q_{L}$", fontsize=14)

This type of identification process is what the LocalDescriptors data outputs may be used for. In the case of Steinhardt OPs, it provides a simple fingerprint for comparing thermalized systems to a known ideal structure to measure their similarity.

For reference, we can also check these values against the Steinhardt class again.

for i, (system, nlist) in enumerate(zip(systems, nlists)):
    steinhardt = freud.order.Steinhardt(l=L)
    steinhardt.compute(system, nlist)
    if np.allclose(steinhardt.particle_order, ld_qls[i][:, L]):
        print("Our manual calculation matches the Steinhardt class!")
Our manual calculation matches the Steinhardt class!
Our manual calculation matches the Steinhardt class!
Our manual calculation matches the Steinhardt class!