Skip to content

Brownian Dynamics

Jump to bottom
Raul edited this page Mar 2, 2021 · 7 revisions

Brownian Dynamics integrators[1]. See BrownianDynamics.cuh

There are several Integrators under the BD namespace, which solve the following differential equation:

dX = K·X·dt + M·F·dt + sqrt(2·kb·T·dt)·B·dW

Where:

  • X - Positions
  • M - Mobility -> M = D/kT
  • K - Shear matrix
  • dW- Brownian noise vector (gaussian numbers with mean=0, std=1)
  • B - B*B^T = M -> i.e Cholesky decomposition B=chol(M) or Square root B=sqrt(M)

The mobility is just a number, so there are no hydrodynamic interactions. See BDHI[5] for solvers with hydrodynamic interactions.

EulerMaruyama

The simplest algorithm, described in [6]. Advances the simulation with the following rule:
X[t+dt] = X[t] + dt(K·X[t]+M·F[t]) + sqrt(2*Tdt)·dW·B
This algorithm has a convergence scaling of 1/2 (dt^0.5).
Given the other alternatives which yield better accuracy with the same cost this solver should be used sparingly.

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:

X[t+dt/2] = X[t] + 0.5dt(K·X[t]+M·F[t]) + sqrt(Tdt)·dW1·B
X[t+dt] = X[t] + dt(K·X[t+dt/2] + M·F[t+dt/2]) + sqrt(Tdt)·B·(dW1 + dW2)

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:

X[t+dt] = X[t] + dt*(K·X[t] + 1.5M·F[t] - 0.5M·F[t-dt]) + sqrt(2Tdt)·B·dW

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):
X[t+dt] = X[t] + dt(K·X[t]+M·F[t]) + sqrt(0.5*Tdt)·B·(dW[t] + dW[t-dt])
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:

  • std::vector<real3> K Shear matrix where K[0] = {Kxx, Kxy, Kxz}, ... Defaults to 0
  • real temperature Temperature of the solvent
  • 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. If particle radius has not been allocated prior construction of EulerMaruyama it will fail.

EXAMPLE

#include"uammd.cuh"
using namespace uammd;
int main(int argc, char *argv[]){
...
//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, pg, sys, par);
...
//Add any interactor[2]
bd->addInteractor(myInteractor);
...
//Take simulation to the next step
bd->forwardTime();
...
return 0;
}

References:
[1] https://github.com/RaulPPelaez/UAMMD/wiki/Integrator
[2] https://github.com/RaulPPelaez/UAMMD/wiki/Interactor
[3] Temporal Integrators for Fluctuating Hydrodynamics. Delong et. al. (2013) Phys. Rev. E 87, 033302.
[4] Brownian dynamics of confined suspensions of active microrollers. Balboa et. al. (2017) J. Chem. Phys. 146; https://doi.org/10.1063/1.4979494
[5] The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics. Leimkuhler et. al. IMA J. Numerical Analysis 36, 1 (2016) https://doi.org/10.1093/imanum/dru056
[6] An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations. Desmond J. Higham. (2001). https://doi.org/10.1137/S0036144500378302

Clone this wiki locally