Three-Body Distribution (g3)

Definition¶

The rooted three-body distribution $g_3$ captures pairwise distance and angular correlations between atomic triplets. For a centre atom at position $\mathbf r_0$ with neighbours at $\mathbf r_1$ and $\mathbf r_2$:

• $r_{01} = \lVert \mathbf r_1 - \mathbf r_0 \rVert$

• $r_{02} = \lVert \mathbf r_2 - \mathbf r_0 \rVert$

• $\phi = \arccos\!\left(\dfrac{(\mathbf r_1 - \mathbf r_0) \cdot (\mathbf r_2 - \mathbf r_0)}{r_{01} \, r_{02}}\right)$

G3Distribution.measure_g3 replicates the supercell up to $r_\text{max}$ and, for every site, accumulates counts into a 4D histogram

over all neighbour pairs $(j, k)$ inside the radial cutoff (with $j \ne k$ when the two are the same species). Radial bins use width $\Delta r = \texttt{r\_step}$; angular bins use $\Delta \phi = \pi / \texttt{phi\_num\_bins}$.

Channel indexing¶

Channels are labelled by species triplets (centre | neigh_1 neigh_2) with species_idx(neigh_1) ≤ species_idx(neigh_2). For a system with $S$ species the channel count is

Examples: SiO₂ ($S = 2$) → 6 channels (Si | Si Si, Si | Si O, Si | O O, O | Si Si, O | Si O, O | O O); SrTiO₃ ($S = 3$) → 18 channels.

g3_index stores the (centre, neigh_1, neigh_2) species indices for each channel, and g3_lookup[c, a, b] maps any ordered species triple to its channel; both orderings of the two neighbours land on the same entry so the caller doesn’t have to sort. The atom-pair $(j, k)$ and $(k, j)$ are added separately when the two neighbour species differ (symmetrising the $(r_{01}, r_{02})$ plane); when they are the same species the diagonal is masked so no atom pairs with itself.

Acceleration backend¶

measure_g3(backend="auto") (the default) uses the numba-parallel kernel. numba is a hard dependency of atomode (installed by pip install atomode), so the fast path is always available; the pure-numpy reference loop remains in the source as a self-check and runs when backend="python" is passed explicitly. Both backends produce bit-identical g3count and g2count (verified by tests/test_g3_numba.py).

The first call pays a one-time JIT compile cost (~1–2 s); subsequent calls hit the cached kernel. The kernel uses prange over origin atoms with per-thread accumulators reduced at the end — no atomic adds, no GIL.

Subsampled measurement (large cells)¶

Measuring g3 directly on a very large cell (100³ Å and beyond) takes minutes even with the numba kernel. measure_g3 accepts a sample_fraction knob that uniformly subsamples origin atoms while keeping the full neighbour set intact — periodicity and histogram shape are preserved while the cost drops linearly:

# ~1% of origins — same statistics as a full measurement at 40³ Å
cell.measure_g3(sample_fraction=0.01, sample_rng_seed=2026)

The subsampled histogram correlates > 0.998 with the full measurement for any sample_fraction ≥ 1/64 (validated on 60³ Å SiO₂), and the bond-peak position is unchanged to three decimals. plot_g2_compare accepts the same kwargs for g(r) overlays on large cells.

Reduced coordinates¶

The random-limit (ideal gas) g3 scales as $r_{01}^2 \, r_{02}^2 \, \sin\phi$. The reduced g3 divides this out:

where $A$ is a per-channel amplitude estimated from the far-field mean. In reduced coordinates $\tilde g_3 \to 1$ in the random limit and crystalline peaks sit as sharp islands above the background. The pairwise byproduct $g_2$ is normalised analogously with the $r^2$ factor divided out.

Coordination shell target¶

CoordinationShellTarget.from_atoms(atoms) extracts first-shell structural targets from a reference crystal using ASE’s periodic neighbor_list at an adaptive cutoff (clamped between 2 Å and $3.8 \times$ the median nearest-neighbour distance). The resulting frozen dataclass exposes (per species pair $(a, b)$):

Triplet channels on the shell target use the same (centre, neigh_1 ≤ neigh_2) indexing (angle_index, angle_lookup) as the g3 object, so downstream consumers can re-use the same channel layout.

Bond restrictions¶

Pairs with zero coordination_target install no bond springs and cannot be bonded during shell_relax. from_atoms zeroes lattice-artefact pairs automatically (auto_filter_lattice_artifacts=True): a species pair whose first-shell peak is not the smallest in either its row or column of pair_peak is a lattice separation through a bridging atom, not a chemical bond. This is what keeps SiO₂ / SrTiO₃ from developing a spurious Si-Si or Ti-Ti bond at the second-shell distance.

Two helpers set the bond graph explicitly instead:

shell_target.with_cross_species_bonds_only()          # zero diagonal: only A–B bonds
shell_target.with_bonded_species_pairs([('Ti', 'O')]) # keep listed pairs, zero the rest

Both return a new CoordinationShellTarget; the original is unchanged.

Angle-spring masking (multi-modal shells)¶

For shells with multiple physically valid bond angles (the 12-coord cuboctahedron at 60°/90°/120°/180°, Sr in SrTiO₃’s SrO₁₂ motif, …) the single angle_mode_deg that gets extracted picks one of the peaks and forcing it strains the others. Two helpers flip entries of angle_enabled_mask to silence specific triplet types without touching their bond-distance springs:

shell_target.with_angle_triplets([('Ti','O','O'), ('O','Ti','Ti')])  # whitelist
shell_target.without_angle_triplets([('Sr','O','O')])                # blacklist

For SrTiO₃ the whitelist form keeps Ti-centered 90° (octahedral) and O | Ti Ti 180° (linear Ti-O-Ti backbone) while silencing every Sr-centered triplet, consistent with Cu FCC’s angle_weight=0 treatment of the cuboctahedron.

Composite shell target¶

CoordinationShellTarget.from_targets({key: target, ...}) stacks several per-chemistry targets into one with a widened species axis. Each input contributes its own virtual-species rows of coordination_target, pair_peak, angle_mode_deg, etc. Cross-target coordination defaults to zero (no bonds form across the virtual boundary), but cross-target repulsion is preserved so different phases don’t overlap. Used for the sp²/sp³ carbon ladder where the two virtual species (sp2_C, sp3_C) share atomic number 6 but carry distinct coordination (3 vs 4) and angle targets (120° vs 109.5°). Per-atom virtual species is assigned at grain-build time via Supercell.generate(..., grain_sources=[...]); see Supercell generation for the grain mechanism.