RSA ensembles across packing fractions

Compare independent monodisperse RSA configurations with Percus–Yevick at three packing fractions. Each analytical curve uses the mean measured RSA fraction; individual runs have the same radius and periodic box size. Shading shows one standard error across independent RSA configurations, not the error of the analytical model or a confidence interval for model discrepancy.

The ensembles are deliberately small for the documentation build. RSA is irreversible and history-dependent, whereas Percus–Yevick is an equilibrium reference; agreement is not expected to be exact. Finite size, radial bins, and the finite ensemble also affect the comparison.

target=0.05; measured range=0.0503–0.0503; independent runs=8
target=0.10; measured range=0.1000–0.1000; independent runs=8
target=0.20; measured range=0.2000–0.2000; independent runs=8

import matplotlib.pyplot as plt
import numpy as np
from PackLab import analytical, monte_carlo, samplers, ureg

radius = 0.2 * ureg.micrometer
domain = monte_carlo.PackingDomain(
    4 * ureg.micrometer,
    4 * ureg.micrometer,
    4 * ureg.micrometer,
    use_periodic_boundaries=True,
)
sampler = samplers.ConstantRadiusSampler(radius, bins=1)
figure, axes = plt.subplots(1, 3, figsize=(12, 3.5), sharey=True)
for axis, target in zip(axes, (0.05, 0.10, 0.20)):
    curves, fractions = [], []
    for seed in range(10, 18):
        options = monte_carlo.RSAOptions()
        options.random_seed = seed
        options.target_packing_fraction = target
        options.maximum_attempts = 100_000
        options.maximum_consecutive_rejections = 20_000
        result = monte_carlo.RSASimulator(domain, sampler, options).run()
        centers, partial_g = result.compute_partial_pair_correlation_function(
            n_bins=60,
            maximum_pairs=100_000,
        )
        curves.append(partial_g[0, 0])
        fractions.append(result.statistics.packing_fraction_geometry)
    mean = np.mean(curves, axis=0)
    standard_error = np.std(curves, axis=0, ddof=1) / np.sqrt(len(curves))
    fraction = np.mean(fractions)
    reference_domain = analytical.PercusYevickDomain(
        size=100 * ureg.micrometer,
        radii=np.array([0.2]) * ureg.micrometer,
        volume_fraction=fraction,
        number_fractions=[1.0],
    )
    reference = analytical.PercusYevickSolver(
        densities=reference_domain.particle_densities_per_radius,
        radii=reference_domain.radii,
        wavenumber="auto",
    ).compute(centers)
    x = centers.to("micrometer").magnitude
    axis.plot(x, mean, label="RSA ensemble mean")
    axis.fill_between(
        x, mean - standard_error, mean + standard_error, alpha=0.25, label="RSA standard error"
    )
    axis.plot(x, reference.g[0, 0], "k--", label="Percus--Yevick")
    axis.set(title=f"mean fraction = {fraction:.3f}", xlabel="separation (µm)")
    axis.grid(alpha=0.2)
    print(
        f"target={target:.2f}; measured range={min(fractions):.4f}–{max(fractions):.4f}; "
        f"independent runs={len(curves)}"
    )
axes[0].set_ylabel("pair correlation g(r)")
axes[0].legend(fontsize="small")
figure.tight_layout()