qdiffusivity.diffusivity

Kernel density estimator for position-dependent diffusivity.

This module provides a kernel-weighted local estimator for the perpendicular (transverse) and parallel diffusivities of nanoconfined molecular dynamics simulations, generalising the per-project script zn-el/analysis/diff_kde/kde_diffusivity.py into a reusable AnalysisBase class.

Theory

The estimator works in u-space, the CDF-uniformised coordinate \(u = F(z) \in [0, 1]\) built from the pooled equilibrium positions. In u-space the equilibrium measure is uniform, so a single global bandwidth is appropriate across the whole confined region (including near the walls, where a z-space bandwidth would over-smooth the sparse adsorption peaks and under-resolve the dense bulk).

  • CDF map. Each starting position \(z\) is mapped to \(u = F(z) \in [0, 1]\) via the empirical CDF (piecewise-linear interpolation through the sorted positions, using the midrank \(p = (i + 0.5)/n\) so the CDF never returns exactly 0 or 1).

  • Perpendicular (transverse) estimator. For each trajectory increment with starting position \(z_j\) and displacement \(\Delta z_j\), the z-space local estimator is \(\hat d_j = (\Delta z_j)^2 / (2\Delta t)\). The kernel localises in u-space (variance equalisation) but the estimator is the z-space MLE, so \(D_\perp\) is obtained directly — no \(\rho^{-2}\) conversion is applied:

    \[\hat{D}_\perp(u_0) = \frac{\sum_j K_h(u_0 - u_j)\,(\Delta z_j)^2/(2\Delta t)} {\sum_j K_h(u_0 - u_j)}.\]

  • Parallel estimator. Only the starting position is mapped to u-space; the parallel displacement \((\Delta x_j, \Delta y_j)\) is not rescaled (parallel motion is unbounded):

    \[\hat{D}_\parallel(u_0) = \frac{\sum_j K_h(u_0 - u_j)\,[(\Delta x_j)^2 + (\Delta y_j)^2] /(4\Delta t)} {\sum_j K_h(u_0 - u_j)},\]

    where \(\langle \Delta x^2 + \Delta y^2\rangle = 4 D_\parallel \Delta t\) gives the \(4\Delta t\) denominator.

  • Kernel. Either the Gaussian kernel \(K_h(x) = \exp(-x^2/(2h^2))/(h\sqrt{2\pi})\) (infinite support, smooth derivatives — the default) or the Epanechnikov kernel \(K_h(x) = \frac{3}{4h}(1 - (x/h)^2)\) for \(|x| < h\) (compact support, no leakage beyond the boundaries).

  • Bandwidth. "auto" uses a Sheather-Jones plug-in (a kernel-appropriate pilot density is built on a fine grid, its second derivative estimated by central differences, and the oracle bandwidth \(h^* = (\|K\|^2 / (N\,\mu_2(K)^2\,\widehat{\int[\hat f'']^2}))^{1/5}\) evaluated); the Silverman rule of thumb is the fallback. "silverman" uses the rule of thumb directly, and a float fixes the bandwidth.

  • Boundary handling. Mirror reflection at \(u = 0\) and \(u = 1\): for each data point at \(u_j\), mirror copies at \(-u_j\) and \(2 - u_j\) are added (with the same \(\hat d_j\)) before kernel evaluation, so kernel mass leaking outside \([0, 1]\) is reflected back.

  • Evaluation grid. n_points cell-centred points uniform in u-space, \(u_m = (m + 0.5)/M\) for \(m = 0, \ldots, M-1\), mapped back to z via \(z_m = F^{-1}(u_m)\).

  • Error bars. The weighted variance of \(\hat d_j\) within the kernel window, divided by the Kish effective sample size \(N_{\mathrm{eff}}(u_0) = (\sum_j K_h)^2 / \sum_j K_h^2\), gives the standard error of the mean at each grid point.

Functions

build_cdf(z_pooled)

Build a smooth CDF from pooled positions.

gaussian_kernel(x, h)

Gaussian kernel \(K_h(x) = \exp(-x^2/(2h^2))/(h\sqrt{2\pi})\).

kde_estimate(u_data, d_data, u_eval, h[, ...])

Kernel-weighted local estimator in u-space with mirror reflection.

select_diff_bandwidth(u_data[, method, kernel])

Select the diffusivity KDE bandwidth.

sheather_jones_bw(u_data[, kernel])

Simplified Sheather-Jones plug-in bandwidth.

silverman_bw(u_data)

Silverman's rule of thumb (Gaussian reference).

Classes

LocalDiffusivityQKDE(atomgroup, *[, dim, ...])

Kernel-weighted local estimator for transverse diffusivity.