freud

freud - Molecular Simulation Analysis

Overview

freud is a Python library for analyzing MD/MC simulation trajectories. It provides parallelized C++ compute classes that accept NumPy arrays or native frame objects and return NumPy arrays.

Core pattern: Create compute object -> call .compute(system, ...) -> access result properties.

Critical: freud boxes are centered at origin. LAMMPS positions must be shifted (see Data Input section).

Data Input

Reading LAMMPS files (use lammpsio, NEVER write a manual parser)

import lammpsio
import freud
import numpy as np

traj = lammpsio.DumpFile("trajectory.lammpstrj")

for snapshot in traj:
    # Convert box (REQUIRED: freud box is centered at origin, LAMMPS is not)
    low, matrix = snapshot.box.to_matrix()
    freud_box = freud.box.Box.from_matrix(matrix)

    # Center positions from LAMMPS origin to freud's centered box
    center = (snapshot.box.low + snapshot.box.high) / 2
    positions = snapshot.position - center

    # Use with any freud compute class
    rdf.compute(system=(freud_box, positions), reset=False)

Wrapped vs unwrapped coordinates: LAMMPS dump columns x y z are wrapped (use for RDF, order params, etc.); xu yu zu are unwrapped (use for MSD). lammpsio reads whichever columns the dump contains into snapshot.position. For pair analyses, ensure positions are wrapped: positions = freud_box.wrap(positions).

For data files: lammpsio.DataFile("init.data").read() returns a Snapshot. To copy topology: lammpsio.DumpFile("traj.lammpstrj", copy_from=data_snapshot).

Reading GSD files (pass frame directly as system)

import gsd.hoomd

traj = gsd.hoomd.open("trajectory.gsd", "r")
frame = traj[-1]  # GSD frames are native system-like objects
rdf.compute(system=frame, reset=False)  # pass directly, no box/position extraction needed

System-like objects accepted by all compute methods

Source	How to pass
GSD frame	Pass directly as system
MDAnalysis Timestep	Pass directly
HOOMD-blue Snapshot	Pass directly
garnett Frame	Pass directly
OVITO DataCollection	Pass directly
lammpsio / raw arrays	(freud.box.Box(...), positions_array) tuple

Box construction

freud.box.Box(Lx, Ly, Lz, xy=0, xz=0, yz=0)  # direct
freud.box.Box.from_box([Lx, Ly, Lz, xy, xz, yz])  # from list/dict/object
freud.box.Box.from_matrix(matrix_3x3)  # from 3x3 matrix
freud.box.Box.from_box_lengths_and_angles(L1, L2, L3, alpha, beta, gamma)
box.wrap(positions)  # wrap positions into box (modifies in-place, also returns)

Quick Reference

density.RDF - Radial Distribution Function

rdf = freud.density.RDF(bins=200, r_max=5.0, r_min=0, normalization_mode="exact")
# normalization_mode: "exact" (default) or "finite_size" (better for small systems)

for frame in trajectory:
    rdf.compute(system=frame, reset=False)  # reset=False accumulates over frames

r = rdf.bin_centers   # (bins,) array
g_r = rdf.rdf          # (bins,) array
n_r = rdf.n_r          # cumulative number density

Tip: Set r_max to at most half the smallest box dimension.

order.Steinhardt - Bond Order Parameters (q_l, w_l)

# Basic q6
ql = freud.order.Steinhardt(l=6)

# Averaged q6 (Lechner-Dellago) - MUCH better for crystal classification
ql_avg = freud.order.Steinhardt(l=6, average=True)

# w_l variant (third-order invariant, more discriminating)
wl = freud.order.Steinhardt(l=6, wl=True, wl_normalize=True)

# Multiple l values at once
ql_multi = freud.order.Steinhardt(l=[4, 6])

ql.compute(system=frame, neighbors=dict(num_neighbors=12))
q_values = ql.particle_order  # (N,) or (N, len(l)) array

Best practice for crystal classification:

• Use average=True (Lechner-Dellago) for FCC/BCC/HCP/liquid discrimination
• Use Voronoi neighbors for Minkowski Structure Metrics: neighbors=voronoi_nlist
• Typical: compute both q4 and q6 with average=True, then scatter plot or threshold

msd.MSD - Mean Squared Displacement

# Window mode (default): better statistics via time-averaging
msd_calc = freud.msd.MSD(mode="window")

# Direct mode: simple displacement from frame 0
msd_calc = freud.msd.MSD(mode="direct")

# For wrapped coordinates, provide box and images for unwrapping
msd_calc = freud.msd.MSD(box=freud_box, mode="window")
msd_calc.compute(positions=pos_array, images=img_array)
# positions: (N_frames, N_particles, 3) array
# images: (N_frames, N_particles, 3) integer array (optional)

msd_values = msd_calc.msd  # (N_frames,) array

# Diffusion coefficient: D = slope / (2 * d) where d = dimensionality
from scipy.stats import linregress
slope, *_ = linregress(time[1:], msd_values[1:])
D = slope / (2 * 3)  # 3D

Important: For unwrapped coordinates (xu, yu, zu in LAMMPS), just pass positions without box/images. For wrapped coordinates, provide box and image flags.

cluster.Cluster - Connected Components

cl = freud.cluster.Cluster()
cl.compute(system=frame, neighbors=dict(r_max=1.5))
labels = cl.cluster_idx  # (N,) cluster label per particle

# Cluster properties
props = freud.cluster.ClusterProperties()
props.compute(system=frame, cluster_idx=cl.cluster_idx)
centers = props.centers       # (n_clusters, 3) center of mass
sizes = props.sizes           # (n_clusters,) particle counts

References

• references/compute_classes.md — Additional compute classes (Hexatic, Nematic, PMFT, GaussianDensity, Interface, StaticStructureFactor, LocalDensity, Voronoi neighbors).

Parallelism Control

import freud

# Context manager (preferred)
with freud.NumThreads(4):
    rdf.compute(system=frame)

# Global setting
freud.set_num_threads(4)
n = freud.get_num_threads()

Common Patterns

Accumulate over trajectory frames

rdf = freud.density.RDF(bins=200, r_max=5)
for frame in trajectory:
    rdf.compute(system=frame, reset=False)  # reset=False is the key
# rdf.rdf is now the frame-averaged result

Per-type RDF (e.g., type A-B correlation)

rdf = freud.density.RDF(bins=200, r_max=5)
for frame in trajectory:
    types = frame.particles.typeid  # or snapshot.typeid with lammpsio
    pos_A = positions[types == 0]
    pos_B = positions[types == 1]
    # query_points=pos_B computes g(r) for B particles around A reference points
    rdf.compute(system=(box, pos_A), query_points=pos_B, reset=False)

Crystal structure classification workflow

# 1. Compute averaged Steinhardt q4 and q6
q4 = freud.order.Steinhardt(l=4, average=True)
q6 = freud.order.Steinhardt(l=6, average=True)
q4.compute(system=frame, neighbors=dict(num_neighbors=12))
q6.compute(system=frame, neighbors=dict(num_neighbors=12))

# 2. Scatter plot q4 vs q6 - clusters identify crystal structures
# Typical reference values (averaged):
#   FCC: q4 ~ 0.19, q6 ~ 0.57
#   BCC: q4 ~ 0.036, q6 ~ 0.51
#   HCP: q4 ~ 0.097, q6 ~ 0.48
#   Liquid: q6 < 0.35

# 3. Or use SolidLiquid for automatic solid/liquid classification
sl = freud.order.SolidLiquid(l=6, q_threshold=0.7, solid_threshold=6)
sl.compute(system=frame, neighbors=dict(num_neighbors=12))
is_solid = sl.cluster_sizes > 0  # boolean per-particle

Common Mistakes

Mistake	Fix
Writing a manual LAMMPS dump parser	Use lammpsio.DumpFile() - it handles all formats and edge cases
Forgetting to center LAMMPS positions	Shift by (box.low + box.high) / 2 before passing to freud
Extracting box/positions from GSD frames	Pass GSD frame directly as system - it's natively supported
Using Steinhardt(l=6) for crystal classification	Use Steinhardt(l=6, average=True) for much better discrimination
Fixed-count neighbors for all analyses	Consider Voronoi neighbors (freud.locality.Voronoi) for parameter-free, adaptive results
reset=True (default) when averaging over frames	Use reset=False to accumulate; the result auto-averages
Using mode="direct" for MSD	Use mode="window" (default) for better statistics via time-window averaging
Setting r_max larger than L/2	RDF/neighbor artifacts occur beyond half the smallest box dimension