Algorithm

Linked-cell k-nearest search in src/knearest.rs. Allen and Tildesley linked cells, a k-heap per source, Chebyshev shells until the k-th neighbour cannot lie outside the visited cube. The search takes no cutoff.

Inputs

• xyz: Cartesian points.
• simbox: periodic parallelepiped (Cell / lc_cell).
• k: neighbours per source. Must be at least 1.
• mask: optional. false / zero drops a point as source and as candidate.
• cell_hint: target bin edge. None or <= 0 uses 3.0, then the value is clamped to the smallest perpendicular box width.

Each source keeps min(k, n_active - 1) neighbours. Packed outputs fill unused slots with -1.

Fold and bin

1. Take fractional coordinates in [0, 1) (Cell::fractional).
2. Store the Cartesian image in the primary cell (Cell::cartesian).
3. Bin in fractional space: floor(s * n*), clamped to [0, n* - 1].
4. Chain each active point onto a linked list per bin.

nx, ny, nz are floor(width / edge), at least 1. cell_min is the smallest of the three actual cell edges (perpendicular width over count). Orthorhombic boxes wrap each lattice direction independently; a sheared box uses Hinv only for this fold.

Shells

For each source, walk integer cell offsets (dx, dy, dz) in Chebyshev shells. Shell reach is the surface max(|dx|, |dy|, |dz|) == reach (reach == 1 also visits the home cell). max_reach is max(nx, ny, nz) / 2 + 1.

Each offset maps to:

• a primary bin, rem_euclid on the cell indices
• a lattice translation, div_euclid counts through Cell::lattice_shift

The pair distance is Cell::dist2_shifted: Cartesian subtract of the folded points plus that shift. The inner loop does not wrap with Hinv.

The same primary bin can appear under more than one wrap. Each wrap is a separate visit. Skipping those repeats (one shift per unique bin) misses images. That construction, and the ortho cheap path, is in MIC and cells.

Heap and stop

A max-heap of size k stores (dist2, index). For k <= 16 it lives on the stack. After each shell, if the heap is full and the worst dist2 is at most (reach * cell_min)^2, no unvisited point can beat the k-th neighbour, and the walk stops.

knearest returns Neighbors rows (indices, dist2), nearest first. knearest_into / lc_knearest write packed indices. knearest_into_d2 / lc_knearest_d2 also write squared distances. knearest_into_many / lc_knearest_many cover a frame-major batch.

knearest_brute is the all-pairs check used by tests and small systems. It calls Cell::dist2 per pair (ortho: three wraps; general: fractional wrap). It is not the production walk.

Parallel

The parallel Cargo feature (on by default) maps sources with rayon. Each source owns its heap.

Device

linkcell::gpu::Workspace is the same walk on a CUDA device. Fold, bin, then a tiled Hillis-Steele exclusive scan (CUB DeviceScan / HOOMD cell offsets). Occupants are stored in cell-major order with an O(1) home slot. The stencil is a precomputed Chebyshev shell table (LAMMPS NStencil, HOOMD d_cell_adj), not a nested 3-D loop. Eight threads share each source and stride occupants; after each shell they merge heaps and apply the host stop. Output is Cabana's 2-D packed n * k list. knearest_into_many covers every frame that shares a cell. Fold uses Hinv; the pair shift is na a + nb b + nc c. k <= 16.