Brownian Dynamics

When the viscous forces are much larger than the inertial forces, i.e. $|\xi\vec{\pvel}| \gg |m\vec{a}|$ , inertial terms becomes irrelevant at very short time scales. Brownian Dynamics (BD)[2]_ takes advantage of such a time scale separation between the particle velocity fluctuations and its displacement and can be interpreted as the overdamped, or non-inertial, limit of Langevin Dynamics. The decorrelation time of the velocity, defined as $\tau_l = m/\xi$ is much faster than the time needed for a particle to move farther than its own size. BD represents the long time limit of the Langevin equation. This is a powerful property of BD, since sampling the probability distributions of the underlying stochastic processes (stemming from the rapid movement of the solute particles) does not require sampling their fast dynamics.

In BD, the coupling between the submerged particles and the solvent is instantaneous. Furthermore, since the particle velocities decorrelate instantly, the only remaining relevant variables are the positions of the colloidal particles.

Neglecting hydrodynamic interactions (see Brownian Hydrodynamics) we will focus on the simple case of a mobility tensor which is non-zero only on the diagonal. In the case of a no-slip sphere of radius $a$ moving inside a fluid with viscosity $\eta$ the bare self-mobility is given by the well-known Stokes drag $M = (6\pi\eta a)^{-1}$ .

The Brownian dynamics equation of motion (without hydrodynamic interactions) are

$d\vec{\ppos} = M\vec{F}dt + \sqrt{2\kT M}d\vec{\noise},$

Where:

$\vec{\ppos}$ - Particle positions ( $\vec{\ppos} = \{\vec{\ppos}_1, \dots, \vec{\ppos}_N\}$ )

$\vec{F}$ - Particle forces

$M = (6\pi \eta a)^{-1}$ - Mobility -> $M = D/\kT$ . Here $\eta$ is the fluid viscosity and $a$ the hydrodynamic radius of the particles.

$d\vec{\noise}$ - Brownian noise vector (gaussian numbers with mean=0, std=1)

Hint

Brownian Dynamics are proper to the so-called Smoluchowski level of description. Particle positions change slowly enough that their velocity decorrelates effectively instantaneously, inertia can be disregarded and positions are the only relevant variable. The dynamics are described by the BD equations of motion (an over-damped Langevin equation). The diffusion time of the particles, $\tau \sim 10^{-3} s$ , dominates this level.

Brownian Dynamics integrators 

There are several Integrators in UAMMD under the BD namespace, which solve the BD equation above:

Note

The code for this module is located in the source Integrator/BrownianDynamics.cuh

EulerMaruyama

The simplest algorithm, described in [6], advances the simulation with the following rule:

$\vec{\ppos}^{n+1} = \vec{\ppos}^n + dt(M\vec{F}^n) + \sqrt{2\kT M dt}d\vec{W}^n$

This algorithm has a convergence scaling of 1/2 ( $dt^{0.5}$ ).

MidPoint

A two step explicit midpoint predictor-corrector scheme (described in [3]). It has a convergence scaling of $dt^4$ at the expense of having twice the cost of a single step method, as it requires to evaluate forces twice per update. Noise has to be remembered as well but in practice it is just regenerated instead of stored. MidPoint updates the simulation with the following rule:

$\vec{\ppos}^{n+1/2} &= \vec{\ppos}^n + \half dt(M \vec{F}^n) + \sqrt{\kT M dt}d\vec{W}^n_1\\ \vec{\ppos}^{n+1} &= \vec{\ppos}^n + dt(M \vec{F}^{n+1/2}) + \sqrt{\kT M dt}(d\vec{W}^n_1 + d\vec{W}^n_2)$

AdamsBashforth

This algorithm, described in [4], uses the forces from the previous step to improve the prediction for the next. It incurs the overhead of storing the previous forces but its computational cost is marginally larger than EulerMaruyama. This algorithm mixes a first order method for the noise with a second order method for the force. It yields better accuracy than EulerMaruyama, although this comes from experience since as of the time of writing no formal work has been done on its weak accuracy. AdamsBashforth updates the simulation with the following rule:

$\vec{\ppos}^{n+1} = \vec{\ppos}^n + dt(1.5M\vec{F}^n - 0.5 M\vec{F}^{n-1}) + \sqrt{2\kT M dt}d\vec{W}^n$

Leimkuhler

Described in [5] (see eq.45). While also a first order method it seems to yield better accuracy than AB and EM. I am not aware of any formal studies of its accuracy. The update rule is very similar to EM but uses noise from two steps (which are generated each time instead of stored):

$\vec{\ppos}^{n+1} = \vec{\ppos}^n + dtM\vec{F}^{n} + \sqrt{2\kT M dt}\half(d\vec{W}^n + d\vec{W}^{n-1})$

Warning

Note that, as stated in [5], while this solver seems to be better than the rest at sampling equilibrium configurations, it does not correctly solves the dynamics of the problem.

Usage

Use it as any other integrator module.

The following parameters are available:

real temperature Temperature of the solvent in units of energy. This is $\kT$ in the formulas.

real viscosity Viscosity of the solvent.

real hydrodynamicRadius Hydrodynamic radius of the particles (same for all particles*)

real dt Time step

bool is2D = false Set to true if the system is 2D

* If this parameter is not provided, the module will try to use the particle’s radius as the hydrodynamic radius of each particle. In the latter case, if particle radii has not been set in ParticleData prior to the construction of the module an error will be thrown.

#include"uammd.cuh"
#include"Integrator/BrownianDynamics.cuh"
using namespace uammd;
int main(){
  //Assume an instance of ParticleData, called "pd", is available
  ...
  //Choose the method
  using BD = BD::EulerMaruyama;
  //using BD = BD::MidPoint;
  //using BD = BD::AdamsBashforth;
  //using BD = BD::Leimkuhler;
  BD::Parameters par;
  par.temperature=1;
  par.viscosity=1;
  par.hydrodynamicRadius=1;
  par.dt=0.01;
  //Optionally you can place a shear matrix, dX = M*F*dt + sqrt(2*D*dt)*dW + K*R
  //par.K = {{1,2,3},{1,2,3},{1,2,3}};
  //Or, if you want to set just one row:
  //par.K[0] = {1,2,3};
  ...
  auto bd = make_shared<BD>(pd, par);
  ...
  //Add any interactor
  bd->addInteractor(myInteractor);
  ...
  //Take simulation to the next step
  bd->forwardTime();
  ...
  return 0;
}

Here, pd is a ParticleData instance.

Note

As usual, any Interactor can be added to this Integrator, as long as it is able to compute forces.

References: