---
name: freud
description: Use when code imports `freud` or `freud-analysis`, or user asks to analyze molecular simulation trajectories (RDF, MSD, order parameters, structure factor, PMFT, Voronoi, clustering). Also use when reading LAMMPS dump/data files for freud analysis via lammpsio.
---

# 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)

```python
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)

```python
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

```python
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

```python
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)

```python
# 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

```python
# 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

```python
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

```python
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

```python
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)

```python
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

```python
# 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 |