Tutorial: Whole-Shape Fit with Bootstrap Uncertainty

This tutorial covers the whole-fit pipeline: an algebraic sphere or cylinder fit to the entire interface shell at once, with optional bootstrap resampling for an angle uncertainty. Pair it with WallDetector.explicit or WallDetector.from_atoms when you have a known wall position from the simulation setup.

---

1. Overview

The whole-fit pipeline differs from the slicing pipeline in two places:

1. The InterfaceExtractor returns a single (N, 3) shell array — every ray’s interface point pooled into one cloud rather than divided into per-slice 2D sub-clouds.

2. The SurfaceFitter.whole() runs one algebraic Taubin fit on that shell — sphere fit for spherical droplets, cylinder fit (algebraic circle in (x, z) with translational invariance along y) for cylindrical droplets. The contact angle follows from the cap geometry cos θ = (z_wall - z_center) / R.

If bootstrap_samples > 0 the fit is repeated on that many bootstrap resamples of the shell, and the resulting standard deviation of the angles is reported as WholeBatchResult.angle_std.

---

2. Wall detector pairing

The full-sphere Fibonacci ray fan (SpaceSampling.rays() with n_rays_sphere=...) emits rays from the droplet COM in all directions, including downward. Those downward rays hit the wall plane and contribute interface points right at the wall, so the lowest shell point lands at z ≈ z_wall. As a result, WallDetector.min_plus_offset() does work with this pipeline.

The two physical alternatives are usually more reliable:

• WallDetector.from_atoms() reads wall atom positions from the trajectory and places the wall plane at the mean of the top atomic layer. Best when the substrate is explicitly modelled.

• WallDetector.explicit() lets you supply z_wall directly — handy when the wall position is known a priori from the simulation setup.

The slicing-mode equivalent of this consideration is much less sensitive because the wall enters only as a horizontal line in each slice’s 2D fit, but for the whole-shape fit the recovered angle depends quite linearly on the wall z, so it pays to be honest about the wall position.

---

3. Example Code

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

filename = "../../tests/trajectories/traj_spherical_drop_4k.lammpstrj"

# --- Step 1: Identify water-oxygen and wall-atom indices ---
wat_find = LammpsDumpWaterFinder(filename, oxygen_type=1, hydrogen_type=2)
oxygen_indices = wat_find.get_water_oxygen_indices(frame_index=0)

# Wall parser: ``liquid_particle_types`` lists what to EXCLUDE
# (i.e. the liquid), leaving the wall atoms.
wall_parser = LammpsDumpWallParser(filename, liquid_particle_types=[1, 2])
carbon_indices = wall_parser.parse(frame_index=0)  # uses internal IDs

# --- Step 2: Build the whole-fit analyzer ---
# Strategies: full-sphere Fibonacci ray fan + whole-fit + from_atoms wall.
analyzer = TrajectoryAnalyzer(
    parser=LammpsDumpParser(filename),
    atom_indices=oxygen_indices,
    droplet_geometry="spherical",
    interface_extractor=InterfaceExtractor(
        sampling=SpaceSampling.rays(
            n_rays_sphere=400,  # 400 rays uniformly over the full sphere
        ),
        density=DensityEstimator.gaussian(density_sigma=3.0),
    ),
    surface_fitter=SurfaceFitter.whole(
        surface_filter_offset=3.0,
        bootstrap_samples=100,  # 100 bootstrap resamples → angle_std
    ),
    wall_detector=WallDetector.min_plus_offset(),
    wall_atom_indices=carbon_indices,  # routed to the wall detector
)

# --- Step 3: Run the analysis ---
results = analyzer.analyze([1])
batch = results.batches[0]
print(f"angle = {batch.angle:.2f}°  ± {batch.angle_std:.2f}° (bootstrap)")
print(f"R     = {batch.popt[3]:.2f} Å, z_wall = {batch.z_wall:.2f} Å")
print(f"rms residual on the shell = {batch.rms_residual:.2f} Å")

---

4. Expected Output

On the water/graphene fixture:

angle = 99.13° ± 0.41° (bootstrap)
R     = 37.54 Å, z_wall = 4.90 Å
rms residual on the shell = 1.27 Å

The recovered R is not the physical droplet radius — it’s the sphere that best fits the recovered shell, which sits slightly inside the actual interface because of the Gaussian-KDE smoothing. The contact angle reading is unaffected by this; the cap geometry is self-consistent with whichever R and z_c the fit produces.

---

5. Comparing whole vs slicing

Both the whole and slicing pipelines on this fixture recover ~95° when paired with WallDetector.min_plus_offset(0)() (the “natural” interface-derived baseline). When you switch the wall detector to from_atoms or explicit at the actual graphene top (≈ 4.9 Å), both pipelines move up to ~99°. The two methods are physically consistent — the slicing pipeline gives you an inter-slice σ for free; the whole pipeline gives you a bootstrap σ instead.

The choice between them is mostly about what kind of uncertainty you want to report:

• per-slice scatter (slicing): tells you whether the droplet is symmetric across the chosen slice axis.

• bootstrap σ (whole): tells you how well-determined the fit is given the shell points; doesn’t expose asymmetry.

---

6. Tips

• Number of rays (n_rays_sphere): 400 is plenty for a ~30 Å droplet; you’d push to 800–1600 for larger droplets where the shell point density per square Å matters.

• Bootstrap samples: 100 is enough to get an angle_std reliable to two significant figures; 1000 will tighten that but costs ~10× more.

• Cylinder droplets still work — pair the whole fitter with SpaceSampling.rays() configured with delta_cylinder and delta_polar instead of n_rays_sphere. The fitter automatically does a 2D circle fit per the cylinder axis convention.

---

7. Alternative configurations

7.1 Cylindrical whole-fit

For a cylindrical droplet, the whole fitter exploits translational symmetry along the cylinder axis: it does one 2D circle fit in the (x, z) plane to the entire shell, treating the \(y\)-coordinate as ignorable. The same recipe, with the cylinder-mode extractor parameters:

analyzer = TrajectoryAnalyzer(
    parser=LammpsDumpParser(cylinder_fixture),
    atom_indices=oxygen_indices,
    droplet_geometry="cylinder_y",  # or "cylinder_x"
    interface_extractor=InterfaceExtractor(
        sampling=SpaceSampling.rays(delta_cylinder=5.0, delta_polar=8.0),
        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),
)

The recovered popt for the cylindrical case is a 4-element array [xc, zc, R, z_wall] rather than the 5-element spherical [xc, yc, zc, R, z_wall] — the cylinder axis \(y\) doesn’t participate in the fit.

7.2 Explicit wall position

When the wall plane is known a priori from the simulation setup (e.g. a 9-3 LJ wall at a specific \(z\)-coordinate, or you simply read the top atomic layer’s \(z\) in an external script):

wall_detector = WallDetector.explicit(z_wall=5.0)

No wall_atom_indices argument needed on the analyzer in this case — the explicit detector ignores any wall-atom data. Useful both for whole-fit and slicing pipelines.

7.3 rays (binning) alternative

Same Fibonacci-sphere geometry, but the density along each ray is estimated with a 1D top-hat histogram instead of a Gaussian KDE:

interface_extractor = InterfaceExtractor(
    sampling=SpaceSampling.rays(n_rays_sphere=400),
    density=DensityEstimator.binning(bin_width=3.0),
)