Generating Very Large Cells¶

Supercell.generate(num_steps=N) runs a full FIRE quench by default, which is the right tool up to ~50³ Å (a few thousand atoms). Beyond that the FIRE step cost grows fast enough that production-scale cells (200³ Å, ~600 k atoms) become impractical without a different strategy.

This page documents the shortcut pipeline Supercell.bond_relax and the matching tools for measuring g(r) and g3 on a large cell.

When to use the shortcut¶

cell side	atoms (SiO₂)	full FIRE	shortcut
40 Å	~5 k	~10 s/regime	not needed
100 Å	~80 k	~2 min	~30 s
200 Å	~600 k	~30 min	~5 min

The shortcut splits the relaxation into two stages: a cheap O(N log N) attract-to-bond-peak + repel-from-hard-core sweep (Supercell.bond_relax) that brings every pair to a physically correct distance, followed by an optional short FIRE finisher to refine bond angles. For grain-based regimes this matches full-FIRE g(r) to within the histogram bin width at ~1/5 the wallclock.

The pipeline¶

import time
import tricor as tc

# Per-regime build.  `kw` is the same regime dict used with
# `Supercell.generate` (grain_size, bond_weight, angle_weight, ...).
is_disordered = kw.get("grain_size") is None

if is_disordered:
    # Liquid / featureless regimes have a broad natural bond
    # distribution that bond_relax cannot reproduce (it pre-pulls
    # bonds to nominal pair_peak and FIRE can't undo that in a short
    # finisher).  Use the full FIRE pipeline for these regimes.
    cell.generate(shell, capture_trajectory=False, show_progress=False, **kw)
else:
    # Ordered regimes: Voronoi tile → bond_relax → short FIRE finisher.
    kw_no_steps = {k: v for k, v in kw.items() if k != "num_steps"}
    cell.generate(shell, num_steps=0,                # Voronoi only
                  capture_trajectory=False, show_progress=False,
                  **kw_no_steps)
    cell.bond_relax(shell, n_iter=40)                # bring pairs to peak
    cell.shell_relax(shell, num_steps=20,            # angle refinement
                     show_progress=False,
                     bond_weight=kw["bond_weight"],
                     angle_weight=kw["angle_weight"],
                     repulsion_weight=kw["repulsion_weight"],
                     hard_core_scale=kw["hard_core_scale"],
                     nonbond_push_scale=kw["nonbond_push_scale"])

Why the regime-aware branch? bond_relax is a hard pull toward the per-pair coordination peak (e.g. Si-O at 1.61 Å). Crystalline + amorphous regimes are happy there. Liquid wants a much looser distribution (peak ~1.78 Å for SiO₂ liquid because of hard_core_scale=1.05 and weak bond springs); pulling the bonds tight first traps FIRE in a basin it can’t escape with 20 finisher steps.

Measuring g(r) and g3 on a big cell¶

Measuring the three-body distribution directly on a 600 k-atom cell would take minutes. measure_g3 accepts a sample_fraction knob that uniformly subsamples origin atoms while leaving the full neighbour set intact — preserves PBC, preserves shape, drops the cost.

# Measure g3 on ~1% of atoms — same statistics as a 40³ Å cell.
cell.measure_g3(sample_fraction=0.01, sample_rng_seed=2026)

The MC sample correlates >0.998 with the full measurement for any sample_fraction ≥ 1/64 (validated at 60³ Å SiO₂). Bond-peak position is identical to three decimals at every sampling rate.

For the inline g(r) overlay viewer, the same kwargs are plumbed through:

gr_cells = {regime: c for regime, c in cells.items()}
tc.plot_g2_compare(gr_cells, r_max=8.0, r_step=0.05,
                   sample_fraction=0.01, sample_rng_seed=2026,
                   title="SiO₂ 200³ Å — full-cell g(r), MC-subsampled")

This measures g(r) on the full periodic cell with subsampled origins — no chunking, no PBC artefacts.

Slicing for the 3D viewer¶

Rendering 600 k atoms inline in Jupyter is slow. Extracting a smaller chunk for the 3D viewer is fine; just make sure the chunk is not treated as itself periodic, or the g(r) measurement on it will wrap atoms across the chunk faces and create spurious sub-Ångström peaks.

from ase.atoms import Atoms
import numpy as np

def chunk_cube(atoms, side=40.0):
    box = np.diag(atoms.cell.array)
    lo = (box - side) / 2
    hi = lo + side
    m = np.all((atoms.positions >= lo) & (atoms.positions < hi), axis=1)
    return Atoms(
        numbers=atoms.numbers[m],
        positions=atoms.positions[m] - lo,
        cell=np.diag([side, side, side]),
        pbc=False,  # critical: the chunk is a SLICE, not a periodic cell
    )

For measuring g(r) on the chunk, the missing PBC truncates the correlation tail past ~side/2. Prefer measuring directly on the full periodic cell with sample_fraction (above).

Optional cleanup: `enforce_hard_core`¶

After FIRE, a few percent of pairs may sit slightly below the hard-core wall (a natural soft tail in the distribution that matches the full-FIRE reference). If you need a strict zero-overlap guarantee for some downstream consumer, call:

cell.enforce_hard_core(shell, n_iter=40)

This is a pure geometric projection — every violating pair gets pushed apart to exactly shell.pair_hard_min[s_i, s_j].

Trade-off: the projection creates a sharp delta in g(r) at the hard wall (atoms pile up at exactly the wall instead of distributing naturally above it). For visual/analysis purposes the natural soft tail is preferable; the projection is opt-in.

Knob reference¶

knob	default	role
`bond_relax(n_iter=...)`	40	cKDTree attract-to-peak + repel-from-hard-core sweeps. 40 is the saturation point for SiO₂.
`shell_relax(num_steps=...)` finisher	20	Short FIRE pass for angle refinement after `bond_relax`. Set to 0 to skip; bond-only convergence is usually already good.
`measure_g3(sample_fraction=...)`	1.0	MC subsampling of origin atoms. Use `~(small_box/big_box)**3` to match a smaller reference cell’s statistics.
`neighbor_update_interval` (liquid)	10	Bump to 60 for liquid (`angle_weight=0`) — the topology rebuild does no useful work when no triplet forces are active.
`enforce_hard_core(n_iter=...)`	40	Opt-in zero-overlap projection. Creates a delta at the wall; usually leave off.