Usage

Hint

See the Available solvers section for a list of available solvers.

Hint

libMobility can compute the thermal drift term, which in general can be approximated using Random Finite Differences (RFD)

\[\boldsymbol{\partial}_{\boldsymbol{X}} \cdot \mathcal{M} = \lim_{\delta \to 0} \frac{1}{\delta} \left\langle \mathcal{M}\left(\boldsymbol{X} + \frac{\delta}{2} \boldsymbol{W} \right) - \mathcal{M}\left(\boldsymbol{X} - \frac{\delta}{2} \boldsymbol{W} \right) \right\rangle_{\boldsymbol{W}}, \hspace{1cm} \boldsymbol{W} \sim \mathcal{N}\left(0,1 \right)\]

The libMobility interface

All solvers present the same set of methods, which are used to set the parameters of the solver, initialize it, and compute the different terms in the equation above.

As an example, the libMobility.SelfMobility solver has the following API:

class libMobility.SelfMobility(self, periodicityX, periodicityY, periodicityZ)

Bases: object

This module ignores hydrodynamic interactions, i.e. the mobility matrix is \(\frac{1}{6\pi\eta a}\mathbb{I}\). This module will only accept open boundary conditions in the three directions.

LangevinVelocities(self, dt, kbt, forces, torques)

Computes the hydrodynamic (deterministic and stochastic) velocities according to the Langevin equation,

\[\begin{split}\boldsymbol{\mathcal{M}}\begin{bmatrix}\boldsymbol{F}\\\boldsymbol{T}\end{bmatrix} + \sqrt{\frac{2k_BT}{\Delta t}}\boldsymbol{\mathcal{M}}^{1/2} \boldsymbol{W} + k_BT \partial_{\boldsymbol{q}}\cdot \boldsymbol{\mathcal{M}}.\end{split}\]

If forces and torques are omitted then only the stochastic part is computed. Calling this function is equivalent to calling Mdot, sqrtMdotW and divM in sequence and applying their respective scalar coefficients, but this can be done more efficiently in a combined fashion in some solvers. By default, this is equvialent to an Euler-Maruyama scheme that uses a random finite difference to compute the divergence term, if necessary, which may not be accurate enough for some applications.

Parameters:

  • dt (float) – Time step \(\Delta t\) for the Langevin equation.

  • kbt (float) – Boltzmann constant times temperature, \(k_B T\), in units of energy.

  • forces (array_like, optional) – Forces acting on the particles.

  • torques (array_like, optional) – Torques acting on the particles. The solver must have been initialized with includeAngular=True.

Returns:

  • array_like – The resulting linear displacements. Returned shape will be the same as the input forces if given, or the positions if no forces are given.

  • array_like – The resulting angular displacements. Returned shape will be the same as the input torques if given, or the forces/positions if no torques are given. This array will be None if the solver was initialized with includeAngular=False.

Mdot(self, forces, torques)

Computes the product of the Mobility matrix with a group of forces and/or torques, \(\boldsymbol{\mathcal{M}}\begin{bmatrix}\boldsymbol{F}\\\boldsymbol{T}\end{bmatrix}\).

It is required that setPositions has been called before calling this function.

At least one of the forces or torques must be provided.

Parameters:

  • forces (array_like, optional) – Forces acting on the particles. Must have shape (N, 3) or (3N), where N is the number of particles.

  • torques (array_like, optional) – Torques acting on the particles. Must have shape (N, 3) or (3N), where N is the number of particles. The solver must have been initialized with includeAngular=True.

Returns:

  • array_like – The linear displacements. The result will have the same format as the forces (or torques if forces are unspecified).

  • array_like – The angular displacements. The result will have the same format as the torques (or forces if torques are unspecified). This will be None if the solver was initialized with includeAngular=False.

clean(self)

Frees any memory allocated by the module.

divM(self, prefactor=1)

Computes the divergence term, \(\boldsymbol{\partial}_\boldsymbol{x}\cdot \boldsymbol{\mathcal{M}}\). It is required that setPositions has been called before calling this function.

Parameters:

prefactor (float, optional) – Prefactor to multiply the result by. Default is 1.0.

Returns:

  • array_like – The resulting linear displacements. Returned shape will be the same as the input forces if given, or the positions if no forces are given.

  • array_like – The resulting angular displacements. Returned shape will be the same as the input torques if given, or the forces/positions if no torques are given. This array will be None if the solver was initialized with includeAngular=False.

initialize(self, viscosity, hydrodynamicRadius, includeAngular = False, tolerance = 0.0001, seed, lanczos_callback)

Initialize the module with a given set of parameters.

Warning

setParameters must be called before this function.

Parameters:

  • viscosity (float) – Viscosity of the fluid.

  • hydrodynamicRadius (float) – Hydrodynamic radius of the particles.

  • includeAngular (bool, optional) – Whether the solver will produce angular velocities. Needed if torques are given. Default is false.

  • tolerance (float, optional) – Tolerance, used for approximate methods and also for Lanczos (default fluctuation computation). Default is 1e-4.

  • lanczos_callback (callable, optional) – Callback function to be called during the Lanczos process. It should take two arguments: the current iteration (int) and the current error (float). Default is None, which means no callback will be used.

setParameters(self, parameter)

Fixes the parameters of the SelfMobility module. For this module, this is just an example on how to set parameters and it is ignored.

Parameters:

parameter (float) – An example parameter that is not used in this module.

setPositions(self, positions)

The module will compute the mobility according to this set of positions.

sqrtMdotW(self, prefactor=1.0)

Computes the stochastic contribution, \({\mathcal{M}}^{1/2} \boldsymbol{W}\), where \(\boldsymbol{\mathcal{M}}\) is the grand mobility matrix and \(\boldsymbol{W}\) is a standard normal Gaussian process.

It is required that setPositions has been called before calling this function.

Parameters:

prefactor (float, optional) – Prefactor to multiply the result by. Default is 1.0.

Returns:

  • array_like – The resulting linear fluctuations. Returned shape will be the same shape as the positions.

  • array_like – The resulting angular fluctuations. Returned shape will be the same shape as the positions.

Example

The following example shows how to use the libMobility.SelfMobility solver to compute the different terms in the equation above.

A libMobility solver is initialized in three steps:

  1. Creation of the solver object, specifying the periodicity of the system.

  2. Setting the parameters proper to the solver.

  3. Initialization of the solver with global the parameters.

"""
Raul P. Pelaez 2021-2025. libMobility python interface usage example.
One of the available solvers is chosen and then, using the same interface the solver is initialized and the mobility is applied.
"""

import numpy as np
import libMobility

# can also import modules individually by name, e.g.
# from libMobility import SelfMobility
# from libMobility import PSE
# ... etc

# Each solver has its own help available, and libMobility itself offers information about the common interface.
# help(libMobility)
# help(libMobility.NBody)
# help(libMobility.PSE)
# help(...


numberParticles = 3
# Each module can be compiled using a different precision, you can know which one with this
precision = np.float32 if libMobility.SelfMobility.precision == "float" else np.float64
pos = np.random.rand(3 * numberParticles).astype(precision)
force = np.ones(3 * numberParticles).astype(precision)

# All solvers present a common interface, but some of them
# need additional initialization via a "setParameters" function
solver = libMobility.SelfMobility(
    periodicityX="open", periodicityY="open", periodicityZ="open"
)

# For the NBody solver, periodicityZ can also be single_wall, e.g.
# solver = libMobility.NBody(periodicityX='open',periodicityY='open',periodicityZ='open')
# solver.setParameters(algorithm="advise")

# solver = libMobility.PSE(periodicityX='periodic',periodicityY='periodic',periodicityZ='periodic')
# solver.setParameters(psi=1, Lx=128, Ly=128, Lz=128,shearStrain=1)

solver.initialize(
    viscosity=1.0,
    hydrodynamicRadius=1.0,
    includeAngular=False,
)
solver.setPositions(pos)

# Some solvers can include angular velocities by changing includeAngular=True in initialize().
# The return is always a tuple of (linear_velocities, angular_velocities).
linear_velocities, _ = solver.Mdot(forces=force)
print(f"{numberParticles} particles located at ( X Y Z ): {pos}")
print("Forces:", force)
print("M*F:", linear_velocities)

linear_fluctuations, _ = solver.sqrtMdotW()
print("M^{1/2} * dW:", linear_fluctuations)

# Some solvers (e.g. SelfMobility) have no thermal drift and return all zeros.
# In general, the thermal drift is non-zero and other solvers (e.g. DPStokes) return a non-zero value.
linear_drift, _ = solver.thermalDrift()
print("Thermal drift:", linear_drift)
solver.clean()