Getting Started

qdiffusivity provides kernel density estimator (KDE) tools for transverse number-density profiles of nanoconfined molecular dynamics simulations, built on MDAnalysis.

Density profiles with the Epanechnikov KDE

The qdiffusivity.TransverseNumDensityQKDE analysis class pools per-frame positions of an AtomGroup along the confined axis and evaluates an Epanechnikov-kernel KDE on a uniform grid spanning the confined region. Kernel mass that would leak beyond the boundaries is folded back by mirror reflection, so the profile is artefact-free at the walls.

import MDAnalysis as mda
from qdiffusivity import TransverseNumDensityQKDE

u = mda.Universe("topology.data", "trajectory.xtc")
ag = u.select_atoms("type 1 2")  # water atoms

kde = TransverseNumDensityQKDE(
    ag,
    dim=2,
    z_bot=10.0,
    z_top=90.0,
    n_points=400,
    grouping="residues",
    bandwidth="auto",
)
kde.run()

# Number density (particles / Å^3):
#   n(z) = (N_total / (n_frames_used * Lx * Ly)) * rho
import numpy as np
Lx, Ly = u.dimensions[:2]
n_density = (kde.n_total /
             (kde.n_frames_used * Lx * Ly)) * kde.rho

The class is an AnalysisBase subclass, so the usual run(start, stop, step) interface applies.

Low-level KDE utilities

The building blocks are also exposed for direct use on pooled position arrays:

from qdiffusivity import (
    epanechnikov_kernel,
    kde_1d,
    select_bandwidth,
)

z_pooled = np.array([...])  # pooled per-frame positions
z_eval = np.linspace(z_bot, z_top, 400)
h = select_bandwidth(z_pooled, z_bot, z_top, method="auto")
rho, n_eff = kde_1d(z_pooled, z_eval, h, z_bot, z_top)

See the API Documentation for full reference.

Diffusivity profiles with the KDE local estimator

The qdiffusivity.LocalDiffusivityQKDE analysis class estimates the perpendicular (transverse) and parallel diffusivities as a function of position along the confined axis. It works in CDF-uniformised u-space, where the equilibrium measure is uniform so a single global bandwidth is appropriate across the whole gap (including near the walls). The perpendicular estimator uses the z-space local estimator \((\Delta z)^2/(2\Delta t)\), kernel-weighted in u-space; the parallel estimator uses \((\Delta x^2+\Delta y^2)/(4\Delta t)\), kernel-weighted by the starting position in u-space. Kernel mass leaking beyond \(u \in [0, 1]\) is folded back by mirror reflection.

import MDAnalysis as mda
from qdiffusivity import LocalDiffusivityQKDE

u = mda.Universe("topology.data", "trajectory.xtc")
ag = u.select_atoms("type 1 2")  # water atoms

kde = LocalDiffusivityQKDE(
    ag,
    dim=2,
    n_points=200,
    bandwidth="auto",
    kernel="gaussian",
)
kde.run()

# D_perp, D_para are in Ų/ps if the trajectory dt is in ps.
# Mask poorly-sampled regions using the Kish effective sample size:
valid = kde.n_eff_perp > 5

Both the Gaussian (infinite support, smooth) and Epanechnikov (compact support, no leakage) kernels are available via kernel="gaussian" or kernel="epanechnikov". The class is an AnalysisBase subclass, so the usual run(start, stop, step) interface applies.

Itô correction

The perpendicular local estimator \((\Delta z)^2/(2\Delta t)\) carries an \(O(\Delta t)\) Itô bias \(\frac{\Delta t}{2}\Phi(z)^2\) where \(\Phi = D(z)\,\rho'(z)/\rho(z)\) in the isothermal (Hänggi–Klimontovich) convention. In wall-bound geometries with adsorption layers this bias is self-suppressing (the \(D^2\) prefactor and the anti-correlation of \(D\) with \(|V'| = |\rho'/\rho|\) make it small — a few % at the walls, negligible in bulk), so it is off by default. To subtract it explicitly, pass ito_correction=True:

kde = LocalDiffusivityQKDE(
    ag, dim=2, n_points=200, ito_correction=True,
)
kde.run()
# kde.ito_bias holds the subtracted (Δt/2) Φ² array; kde.D_perp
# is the bias-corrected perpendicular diffusivity (clipped >= 0).

The parallel estimator has zero Itô bias (no parallel drift) and is unaffected by this option.

Binned (histogram-style) profiles

For users who prefer histogram-style profiles over kernel smoothing, qdiffusivity.binned provides CDF-binned counterparts to the KDE classes. Binning is always in u-space (CDF-uniformised), so bins are naturally finer where the particle density is high and every bin receives a comparable number of samples — the same equal-population strategy as the project’s quantile scripts.

The bins parameter accepts:

  • int — N uniform u-space bins with cloud-in-cell (CIC) assignment (each sample is linearly split between the two nearest bin centres, avoiding bin-edge discontinuities).

  • “quantile” — shortcut for 30 uniform u-space bins (CIC).

  • array_like — explicit u-space edges in [0, 1] with hard assignment (standard histogram behaviour).

from qdiffusivity import (
    TransverseNumDensityQBinned,
    LocalDiffusivityQBinned,
)

# Density profile, 30 quantile bins (CIC):
binned_dens = TransverseNumDensityQBinned(
    ag, dim=2, z_bot=10.0, z_top=90.0, bins="quantile",
)
binned_dens.run()

# Diffusivity profile, 20 bins, with Ito correction:
binned_diff = LocalDiffusivityQBinned(
    ag, dim=2, bins=20, ito_correction=True,
)
binned_diff.run()

Both classes are AnalysisBase subclasses, so the usual run(start, stop, step) interface applies. The diffusivity class supports the same ito_correction keyword as the KDE version.

Two-pass parallelization

The diffusivity classes accept a density_result keyword — a pre-computed TransverseNumDensityQKDE (or TransverseNumDensityQBinned) whose CDF closures are reused for the u-space mapping. If None (default), a density analysis is run internally. Passing a pre-computed result enables a two-pass parallelization strategy:

# Pass 1: density profile (parallelizable via split-apply-combine).
dens = TransverseNumDensityQKDE(ag, dim=2, n_points=200)
dens.run()

# Pass 2: diffusivity using the pre-computed CDF.
# Each frame's u-mapping is stateless, so this pass is also
# parallelizable (the minimum-image convention replaces NoJump).
kde = LocalDiffusivityQKDE(
    ag, dim=2, n_points=200, density_result=dens,
)
kde.run()

All displacements across periodic boundaries use the minimum-image convention (stateless, frame-local), not the stateful NoJump transformation, so each frame can be processed independently.