libMobility: GPU Hydrodynamic Simulations

• libMobility is an open-source Python library that simulates hydrodynamic interactions at the Smoluchowski level using the RPY framework.
• It employs a unified Python interface over CUDA-enabled solvers to efficiently compute deterministic mobility, stochastic Brownian increments, and thermal drift corrections.
• The library supports diverse geometries—from free-space and periodic to confined systems—enabling large-scale, accurate simulations in colloidal and microfluidic applications.

libMobility is an open-source, GPU-accelerated Python library for simulating hydrodynamic interactions in particulate systems at the Smoluchowski level, that is, at the level of overdamped Langevin or Brownian dynamics in which the fluid degrees of freedom have been eliminated and particle motion is governed by a mobility operator. It is designed around the Rotne-Prager-Yamakawa (RPY) level of hydrodynamics and provides CUDA-enabled solvers for deterministic displacements, stochastic Brownian increments, and thermal-drift corrections across several geometries and boundary conditions. The library’s central contribution is a unified Python interface over geometry-specific GPU solvers, so that large-scale hydrodynamic simulation need not be rewritten separately for free space, periodic domains, or confined wall-bounded settings (Fish et al., 2 Oct 2025).

1. Scientific scope and modeling level

At the scientific level, libMobility addresses a standard difficulty of micron-scale fluid-particle modeling: hydrodynamic interactions are long-ranged, many-body, geometry-dependent, and frequently stochastic. In the regime targeted by the library, the fluid is not evolved as an independent dynamical system; instead, it is integrated out, leaving a particle-level description in which forces and torques are mapped to translational and rotational motion by a grand mobility operator. This is the Smoluchowski or overdamped limit, so neither fluid inertia nor particle inertia is resolved (Fish et al., 2 Oct 2025).

The intended physical objects are radially symmetric spherical particles. This assumption is technically consequential. Because the particles are spheres, the mobility is independent of orientation, the orientation divergence vanishes, and only the position-space divergence can contribute to thermal drift. The framework therefore targets Brownian suspensions in which translational motion, rotational motion, translation-rotation coupling, and detailed-balance-preserving stochastic forcing are all relevant, but shape-anisotropic particle hydrodynamics are outside its stated scope.

The library is positioned for colloidal suspensions, biophysics, microbiology, microfluidics, polymeric fluids, and materials science. In this context, “mobility” has its classical Stokesian meaning: the operator relating applied forces and torques to overdamped particle velocities. A common misconception is to read the name as referring to human mobility or network mobility; here it denotes hydrodynamic mobility at zero Reynolds number.

2. Continuum formulation and Brownian dynamics

The underlying continuum description is fluctuating Stokes flow for a zero-Reynolds-number fluid. The governing equations are written as

$$\eta\nabla^2\mathbf{v} = \nabla p - \mathbf{f} - \nabla\cdot\mathcal{Z}, \qquad \nabla\cdot\mathbf{v} = 0,$$

where $\mathbf{v}$ is the fluid velocity, $p$ is pressure, $\eta$ is viscosity, $\mathbf{f}$ is the force density, and $\mathcal{Z}$ is a fluctuating stress. The random stress has zero mean and covariance

$$\langle \mathcal{Z}_{ik}(\mathbf{x},t)\mathcal{Z}_{jm}(\mathbf{x}',t')\rangle =2k_B T\eta(\delta_{ij}\delta_{km}+\delta_{im}\delta_{kj})\delta(\mathbf{x}-\mathbf{x}')\delta(t-t'),$$

so that the fluctuation-dissipation theorem is satisfied (Fish et al., 2 Oct 2025).

Particle forces $\mathbf{F}_i$ and torques $\mathbf{T}_i$ enter through regularized immersed-boundary-style couplings, and particle translational and angular velocities are obtained from corresponding kinematic constraints. Formal elimination of the fluid through the Stokes Green’s function yields the grand mobility matrix $\mathbf{M}$, a $\mathbf{v}$0 operator when both translations and rotations are included. In compact deterministic form,

$\mathbf{v}$1

with stacked translational and angular velocities in $\mathbf{v}$2 and stacked forces and torques in $\mathbf{v}$3.

The Brownian dynamics equation used by the library is

$\mathbf{v}$4

Its three contributions are the deterministic mobility drift, the stochastic increment, and the thermal drift term required for detailed balance when the mobility depends on configuration. The covariance of the Brownian increment is

$\mathbf{v}$5

A central computational implication follows immediately. Direct assembly and factorization of the full mobility matrix is generally infeasible at large $\mathbf{v}$6; with translations and rotations, the matrix dimension is $\mathbf{v}$7, and a naive square-root computation costs $\mathbf{v}$8. libMobility is therefore organized around applying $\mathbf{v}$9, sampling $p$0, and evaluating $p$1 without constructing dense matrix factorizations (Fish et al., 2 Oct 2025).

3. Solver modules and numerical methods

The library exposes several solver backends through a common abstraction. The design principle is not “build the matrix,” but “apply the operator.” This makes geometry-specific kernels interchangeable at the Python level while preserving highly specialized CUDA implementations underneath.

Solver	Supported setting	Asymptotic cost
SelfMobility	Fully open boundaries; neglects interparticle hydrodynamic interactions	$p$2
PSE	Triply periodic bulk systems	$p$3 in particle count
NBody	Free-space RPY or single-wall systems	$p$4
DPStokes	Doubly periodic systems with $p$5 open, single_wall, or two_walls	$p$6

SelfMobility is the simplest backend. It neglects interparticle hydrodynamic interactions and applies single-particle Stokes drag independently, with translational self-mobility $p$7. It is intended as a verification baseline or for dilute settings in which coupling is negligible.

For triply periodic bulk suspensions, libMobility includes the GPU implementation of the Positively Split Ewald algorithm of Fiore et al. PSE is the preferred bulk workhorse. It uses a real-space/reciprocal-space decomposition of the Stokes Green’s function: the real-space part is evaluated in a purely Lagrangian particle-neighborhood form, whereas the reciprocal-space part is computed on a grid. PSE is therefore explicitly a hybrid Eulerian-Lagrangian method. The paper states that it scales nearly linearly with system size and exactly linearly with number of particles, while stochastic forcing is sampled directly in Fourier space at the cost of “one forward and one inverse FFT” (Fish et al., 2 Oct 2025).

The NBody backend is a direct Green’s-function-based brute-force evaluation of pairwise RPY mobilities. Its $p$8 complexity makes it most appropriate for relatively small systems, stated as less than about $p$9K particles, or for very large physical domains where grid-based methods become inefficient. In the $\eta$0 direction it supports either open, giving free-space RPY, or single_wall, giving the Swan-Brady wall-corrected kernel. A numerical subtlety arises because the Swan-Brady kernel is not symmetric positive definite when particles overlap the wall. Since positive definiteness is critical for Brownian noise generation, libMobility incorporates the damping-matrix approach described in Appendix A of the authors’ earlier magnetic-roller work, smoothly damping the mobility to zero near the wall while preserving SPD.

For doubly periodic confined systems, the DPStokes backend implements the spectral method of Hashemi et al. It is asymptotically linear and is intended for dense monolayers and thin-film or slit geometries, requiring periodic $\eta$1 while allowing $\eta$2 to be open, single_wall, or two_walls. In practice it is the confined analogue of PSE (Fish et al., 2 Oct 2025).

Two generic algorithmic layers sit above the deterministic mobility action. For stochastic increments, if a solver lacks its own fluctuation machinery, libMobility uses the Lanczos algorithm to compute $\eta$3 through repeated applications of $\eta$4 without forming $\eta$5 explicitly; the paper states that Lanczos is currently used for DPStokes and NBody. For thermal drift, when specialized code is unavailable, the library uses Random Finite Differences,

$\eta$6

again reducing the problem to repeated mobility applications.

4. Geometries, API, and software architecture

The supported domain descriptors in the constructor are periodic, open, single_wall, and two_walls, realized differently by different backends. The presently supported physical settings are fully open or free-space systems, triply periodic bulk systems, single-wall systems, and doubly periodic slit-like systems with one or two walls. This geometrical distinction is not superficial. In translationally invariant environments such as free space or fully periodic bulk, thermal drift often vanishes; in confined geometries such as bottom-wall or slit-channel systems, the mobility depends on height and the drift term must generally be included (Fish et al., 2 Oct 2025).

The Python API exposes the three Brownian-dynamics primitives separately:

• Mdot: deterministic mobility action
• sqrtMdotW: stochastic displacement
• divM: thermal drift

These appear alongside constructor and state-management methods such as setParameters(...), initialize(...), and setPositions(pos). The interface returns both translational and angular outputs, and the includeAngular flag allows a translational-only mode when rotational blocks are unnecessary.

A representative usage pattern given in the paper is:

$\mathbf{f}$3

Architecturally, libMobility is a Python front end over CUDA-enabled backend solvers and draws on GPU-optimized numerical modules from libraries such as UAMMD. Interoperability is handled through DLPack: any array object supporting the buffer protocol can be passed, including NumPy arrays, CuPy arrays, PyTorch tensors, and JAX tensors. If the input is on the CPU, it is copied to the GPU automatically; if it is already resident in GPU memory, transfers are avoided. The output type generally matches the input type. This design makes the package usable inside standard Python scientific workflows while keeping the heavy operator applications on the GPU (Fish et al., 2 Oct 2025).

5. Performance, validation, and scientific use

libMobility is primarily a GPU-focused library and currently requires an NVIDIA GPU. Prebuilt packages support Linux, and source builds on Windows are reported as possible. CUDA acceleration is used for deterministic mobility application, FFT-based reciprocal-space calculations, pairwise mobility evaluation, iterative square-root sampling, and iterative drift estimation. Benchmarks for Mdot are reported on an NVIDIA RTX A6000. The qualitative scaling trends match the asymptotic analysis: PSE and DPStokes show near-linear or asymptotically linear behavior, NBody follows the expected $\eta$7 trend, the wall-corrected NBody kernel is somewhat more expensive than free-space NBody, and small systems are dominated by GPU overhead and underutilization (Fish et al., 2 Oct 2025).

The library is validated in three categories. Deterministic mobility validation compares computed displacements against analytical or benchmark results, including NBody pair mobilities against exact RPY and Swan-Brady kernels, and PSE and DPStokes against literature spectral or Ewald results. Fluctuation-dissipation validation tests whether sampled fluctuations reproduce the covariance implied by $\eta$8; for this purpose, the full mobility matrix is explicitly factorized by SVD in small validation problems, and agreement is assessed statistically with a Kolmogorov-Smirnov test. Thermal drift validation compares one RFD implementation against an alternate RFD-based estimate and checks agreement on average over many realizations (Fish et al., 2 Oct 2025).

The scientific examples are broader than a minimal API demonstration. In wall-bounded active matter, the NBody backend is used for torque-driven rollers near a wall, including deterministic and stochastic dynamics for $\eta$9 particles, reproducing the formation of “critters” and the delayed instability onset in the stochastic case. In passive-colloid simulations above a wall and in free space, the library is used to extract self-diffusion and intermediate scattering functions; the study reports that above a wall the collective diffusion coefficient appears to plateau at small $\mathbf{f}$0, unlike the known $\mathbf{f}$1 divergence in quasi-2D free space. In rheology, the PSE mobility operator is coupled to a rigid multiblob method and a GMRES resistance solve, with the largest periodic multiblob simulations using roughly half a million blobs, solving in 73 iterations, and running at about 9 iterations per second. In electro-osmotic microfluidics, DPStokes is combined with immersed-boundary boundary blobs and an electrostatic slip model in a slit channel, producing flow fields and velocity profiles that match experiment closely (Fish et al., 2 Oct 2025).

These examples clarify the intended scientific role of the software. libMobility is not merely an implementation of isolated kernels; it is a reusable mobility engine for overdamped particulate hydrodynamics, including use inside larger multiblob, rheological, and multiphysics workflows.

6. Limitations, availability, and significance

The package has clear domain restrictions. It operates at the Smoluchowski or overdamped level and does not resolve fluid or particle inertia. The theory and solvers described are for spherical particles at the RPY level. Hydrodynamics are mainly far-field regularized; near-field lubrication is not automatically included in the core solvers, and the passive-colloid example adds lubrication corrections externally. The solvers do not include stresslets, which limits direct rheology calculations unless the library is coupled to another method such as rigid multiblobs. Some stochastic machinery depends on $\mathbf{f}$2 being symmetric positive definite, geometry support is solver-dependent, and the present implementation is NVIDIA/CUDA dependent rather than portable across accelerator backends (Fish et al., 2 Oct 2025).

Availability is straightforward. The repository is https://github.com/stochasticHydroTools/libMobility, the documentation is https://libmobility.readthedocs.io, the license is MIT, and installation is available through conda install -c conda-forge libmobility. These distribution details reinforce the paper’s presentation of libMobility as production-oriented scientific software rather than a one-off code release.

The main contribution of libMobility is therefore not a single new hydrodynamic kernel but a unified, modular software layer over several specialized hydrodynamic solvers at the Smoluchowski/RPY level. Its significance lies in the combination of three elements: a common abstraction for deterministic mobility, Brownian noise, and thermal drift; production-ready CUDA solvers across open, periodic, and confined geometries; and generic stochastic and drift algorithms that convert any suitable SPD mobility action into a Brownian hydrodynamics solver. A plausible implication is that the library lowers the barrier between method development and large-scale application: geometry-specific, GPU-intensive hydrodynamic operators become accessible through one consistent Python interface while preserving the technical distinctions that matter physically and numerically.