Tutorial: Grid-Based Space Sampling for Interface Extraction

This tutorial covers the grid-based space sampling for interface extraction — SpaceSampling.grid(). It is an alternative to the ray-fan extractors used in the Tutorial: Contact Angle Analysis (Slicing Pipeline) and Tutorial: Whole-Shape Fit with Bootstrap Uncertainty: instead of locating the interface as the half-density point of a 1D tanh fit along each ray, it evaluates a density at each cell of a fixed grid and recover the interface as the iso-density contour at the half-bulk level.

In slicing mode the grid sampling iterates per slice — per azimuthal angle for spherical droplets, per axial step for cylinder droplets — so the downstream SurfaceFitter.slicing sees one contour per slice and can expose per-slice asymmetry, exactly like the rays variants.

1. When to pick grid over rays?

Both extractors plug into the same TrajectoryAnalyzer and produce the same downstream result objects, so the choice is mostly about how the noise/cost trade-off lands on your system:

  • Ray fans sample density along a small number of well-chosen directions; each ray’s 1D tanh fit is cheap. Best when atom statistics per frame are high and you want sub-frame resolution.

  • Grids estimate density on every cell of a fixed mesh, then trace an iso-contour. Closer to the “average over many frames” intuition; the per-cell density gets smoother as more frames are pooled.

The grid sampling requires scikit-image for the iso-contour tracing (marching squares in 2D, marching cubes in 3D). Install via the grid3d extra:

pip install wetting-angle-kit[grid3d]

2. Worked example: grid sampling with gaussian density + slicing fit

A spherical droplet, with per-azimuthal-slice 2D density grids in the (s, z) plane — same density estimator as rays (Gaussian), just sampled on a fixed grid rather than along rays:

from wetting_angle_kit.analysis import (
    DensityEstimator,
    InterfaceExtractor,
    SpaceSampling,
    SurfaceFitter,
    TrajectoryAnalyzer,
    WallDetector,
)
from wetting_angle_kit.analysis.temporal import TemporalAggregator
from wetting_angle_kit.parsers import LammpsDumpParser, LammpsDumpWaterFinder

filename = "../../tests/trajectories/traj_spherical_drop_4k.lammpstrj"
oxygen_indices = LammpsDumpWaterFinder(
    filename, oxygen_type=1, hydrogen_type=2
).get_water_oxygen_indices(frame_index=0)

# 2D grid for each slice plane: ``s`` (in-plane radial) spans
# ``[xi_0, xi_f]`` symmetrically around the slice centre to cover
# the full diameter; ``z`` stays in the lab frame.
grid_params = {
    "xi_0": -40.0,
    "xi_f": 40.0,
    "dx": 3.0,  # 3 Å cells in s
    "zi_0": 0.0,
    "zi_f": 40.0,
    "dz": 1.6,  # 1.6 Å cells in z
}

analyzer = TrajectoryAnalyzer(
    parser=LammpsDumpParser(filename),
    atom_indices=oxygen_indices,
    droplet_geometry="spherical",
    interface_extractor=InterfaceExtractor(
        sampling=SpaceSampling.grid(
            grid_params=grid_params,
            delta_azimuthal=20.0,  # 9 azimuthal slices
        ),
        density=DensityEstimator.gaussian(density_sigma=2.0),
    ),
    surface_fitter=SurfaceFitter.slicing(surface_filter_offset=3.0),
    wall_detector=WallDetector.min_plus_offset(offset=0.0),
    temporal_aggregator=TemporalAggregator(batch_size=1),
)
batch = analyzer.analyze([1]).batches[0]
print(
    f"Angle (grid (Gaussian) + slicing): {batch.angle:.2f}° "
    f"± {batch.angle_std:.2f}° across {len(batch.per_slice_angles)} slices"
)

3. Histogram alternative: grid sampling with binning density

Same per-slice iteration, but the density estimator is a top-hat histogram of atoms within the slab |perp| ≤ dx / 2 of the slice plane. Numerically cheaper than the KDE; intrinsically noisier because only atoms in the slab contribute (not all atoms along the slice direction the way they do for rays (binning)). Use coarser cells (thicker slab) than for grid (Gaussian):

from wetting_angle_kit.analysis import (
    DensityEstimator,
    InterfaceExtractor,
    SpaceSampling,
)

grid_params = {
    "xi_0": -40.0,
    "xi_f": 40.0,
    "dx": 8.0,  # thick slab
    "zi_0": 0.0,
    "zi_f": 40.0,
    "dz": 3.0,
}
extractor = InterfaceExtractor(
    sampling=SpaceSampling.grid(
        grid_params=grid_params,
        delta_azimuthal=60.0,  # fewer slices → more atoms per slab
    ),
    density=DensityEstimator.binning(),
)

The slab thickness perpendicular to each slice plane is dx, so refining the in-plane grid also thins the slab. For systems with limited atom statistics per slab, the answer is either coarser cells or fewer slices, not a finer grid.

4. 3D iso-surface for the whole-fit

The grid sampling also works in whole-fit mode for spherical droplets — the 2D density grid is replaced by a 3D one, and the half-bulk iso-surface is traced via marching cubes. Whole mode takes no delta_azimuthal / delta_cylinder:

grid_params_3d = {
    "xi_0": -30.0,
    "xi_f": 30.0,
    "dx": 2.5,
    "yi_0": -30.0,
    "yi_f": 30.0,
    "dy": 2.5,
    "zi_0": 0.0,
    "zi_f": 35.0,
    "dz": 2.0,
}

analyzer = TrajectoryAnalyzer(
    parser=LammpsDumpParser(filename),
    atom_indices=oxygen_indices,
    droplet_geometry="spherical",
    interface_extractor=InterfaceExtractor(
        sampling=SpaceSampling.grid(grid_params=grid_params_3d),
        density=DensityEstimator.gaussian(density_sigma=3.0),
    ),
    surface_fitter=SurfaceFitter.whole(
        surface_filter_offset=3.0,
        bootstrap_samples=100,
    ),
    wall_detector=WallDetector.min_plus_offset(offset=0.0),
)
batch = analyzer.analyze([1]).batches[0]
print(
    f"Angle (grid (Gaussian) + whole-fit): "
    f"{batch.angle:.2f}° ± {batch.angle_std:.2f}°"
)

Three notes on the 3D case:

  • xi/yi are in the droplet-centred frame (the per-frame COM is subtracted before binning); zi stays in the lab frame so the wall position keeps its physical meaning.

  • grid + whole-fit works for both spherical and cylinder droplets. For a cylinder, the user must pick yi_0 / yi_f to span the full cylinder axis (typically [-box_y / 2, +box_y / 2]); the centred-grid convention puts the droplet COM at the midpoint along y, so a symmetric range covers the whole ridge. The fitter projects the 3D shell onto the (x, z) plane and does a 2D circle fit by translational invariance along y.

  • Marching cubes can be slow on dense 3D grids; if performance matters, start with 2–3 Å cells and only refine if the recovered angle is grid-resolution-limited.

5. Tips

  • Grid bounds: pick xi_f, yi_f, zi_f so the full droplet fits comfortably inside the grid (signed xi_0 for the slicing case so the slice spans the full diameter). The iso-contour tracer can’t extrapolate.

  • Cell sizes: dx controls in-plane horizontal resolution; dz controls vertical. The range bounds are honoured exactly and the cell width is rounded to fit, so the effective cell size may differ slightly from the value you pass.

  • Comparison plot: run the same trajectory through both rays (Gaussian) and grid (Gaussian) and check the two angles agree within method-dependent tolerance (a few degrees on 4k-atom droplets). If they diverge by more than ~8°, one of them is misconfigured (most often the grid bounds are too tight or surface_filter_offset is too small).

  • grid + binning slab thickness: the slab perpendicular to each slice equals dx. If you see a noisy iso-contour, thicken it (larger dx) before reaching for grid (Gaussian).