Langevin Dynamics

At the Langevin level of description, the dynamics of the solvent molecules (e.g. water) are orders of magnitude faster than the solute particles, leaving the positions and momenta of the solute particles as relevant variables. The dynamics of the solute particles are governed by a so-called Langevin equation, a manifestation of the Fokker-Plank equation that is often described as “Molecular Dynamics with a thermostat”.

We can write the Langevin stochastic differential equation as

$m\, d\vec{\pvel} = \underbrace{\vec{F}}_{\vec{\partial}_{\vec{\ppos}}U(\ppos)}dt - \overbrace{\xi\vec{\pvel}}^{\text{Drag}}dt + \underbrace{\vec{\beta}}_{\text{Fluctuations}},$

which can be interpreted as a form of Molecular Dynamics coupled with a thermostat so that noise and drag forces are balanced according to the fluctuation-dissipation relation, guaranteeing the correct thermalization of the system (NVT ensemble). Here $\vec{\pvel} = \dot{\vec{\ppos}}$ are the velocities of the particles.

On the other hand, the LD equation is an expression of the equivalent stochastic differential equation (SDE) for the Fokker-Plank equation with the particle momenta and positions as the relevant variables. The first term in the Langevin SDE , $\vec{F}$ , is the sum of the conservative forces acting on particle $i$ (the ones coming from the underlying potential energy landscape), these interactions can be steric, electrostatic, bonded, etc. In UAMMD, forces are provided to the Integrator modules via Interactors. The second term represents the drag exerted by the solvent particles in the form of a dissipative force. Finally, the third term represents the fluctuations produced by the fast and constant collisions of the solvent particles. Here $\vec{\beta}$ is a random increment which must be in fluctuation-dissipation balance with the friction term.

The friction coefficient, $\xi$ , is related with the damping rate as $\gamma = \xi/m$ , which represents the decorrelation time of the velocity, $\tau = m/\xi$ . Additionally, $\xi$ can be formally derived from a Green-Kubo relation involving the integral of the solvent-solute force time-correlation. Its value is often approximated by the Stokes (macroscopic) value $\xi=6\pi\eta a$ , of a particle of radius $a$ in a solvent with viscosity $\eta$ .

The Fokker-Planck formalism teaches us that when a system dissipates energy, for instance, due to a force like the drag force in the Langevin SDE, the fluctuation-dissipation theorem states that there must be an opposite process that reintroduces this energy via thermal fluctuations.

In this case, fluctuation-dissipation enforces a fluctuating term with zero mean,

$\left\langle\beta\right\rangle = 0,$

which is uncorrelated in time (since we are assuming the time interval is small enough to consider the transport coefficients as constants) with standard deviation given by

$\left\langle\beta(t)\beta(t')\right\rangle = 2\xi\kT dt\delta(t-t').$

We can then write the Langevin equation as

$m\, d\vec{\pvel} = \vec{F}dt - \xi\vec{\pvel}dt + \sqrt{2\xi\kT}\vec{\noise},$

Where we have defined $\vec{\beta} := \sqrt{2\xi\kT}\vec{\noise}$ .

Technically, the restrictions for $\vec{\beta}$ would be compatible with any random distribution for the noise, $\vec{\noise}$ , that has zero mean and variance $dt$ . However, as the noise term comes from the action of countless independent random variables (collisions with the solvent particles) the central limit theorem applies. Thus, it is convenient to choose $\vec{\noise}$ as a Gaussian white noise (i.e. a Wiener process).

Hint

At the Langevin level of description, particles move so slowly compared to the solvent characteristic times that hydrodynamic interactions are effectively instantaneous. Only the positions and momenta of the particles remain as relevant variables (Lagrangian description once again), which obey a Langevin equation. The characteristic decorrelation time of the particle’s velocities governs this scale with $\tau\sim 10^{-5}s$ .

UAMMD exposes a Verlet-like algorithm to numerically integrate the Langevin SDE called Grønbech-Jensen [1]

The update rule, taking the simulation from the time step $n$ (time $t=\dt n$ ) to the next, $n+1$ , can be summarized as

$\vec{\ppos}^{n+1} &= \vec{\ppos}^n + b \dt \vec{\pvel}^n + \frac{b\dt^2}{2m}\vec{F}^n + \frac{b\dt}{2m}\vec{\beta}^{n+1}\\ \vec{\pvel}^{n+1} &= a\vec{\pvel}^n + \frac{\dt}{2m}\left(a\vec{F}^n + \vec{F} ^{n+1}\right) + \frac{b}{m}\vec{\beta}^{n+1}$

Where

$b:&=\frac{1}{1+\frac{\xi\dt}{2}}\\ a:&=b \left(1-\frac{\xi\dt}{2}\right)$

Note that, since we need the force at $t+\dt$ to compute the velocities at $t+\dt$ the forces have to be a function of the positions only (and not velocities).

Note

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

Usage

Grønbech-Jensen [1] is available as an Integrator.

The following parameters are available:

real temperature Temperature of the solvent in units of energy. This is $\kT$ in the formulas.
real friction Friction, $\xi$ .
real dt Time step.
real mass = -1 Mass of all the particles. If >0 all particles will have this mass, otherwise the mass for each particle in ParticleData will be used. If masses have not been set in ParticleData the default mass is 1 for all particles.
bool initVelocities=true Modify starting velocities to ensure the target temperature from the start. When false the velocities of the particles are left untouched at initialization. The default is true and sets particle velocities following the botzmann distribution.
bool is2D = false Set to true if the system is 2D

#include"uammd.cuh"
#include"Integrator/VerletNVT.cuh"
using namespace uammd;
int main(){
  //Assume an instance of ParticleData, called "pd", is available
  ...
  using NVT = VerletNVT::GronbechJensen;
  NVT::Parameters params;
  params.temperature = 1.0;
  params.dt = 0.1;
  params.friction = 1.0;
  auto verlet = std::make_shared<NVT>(pd, params);
  ...
  //Add any interactor
  verlet->addInteractor(myInteractor);
  ...
  //Take simulation to the next step
  verlet->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.

Dissipative Particle Dynamics

One of the most popular techniques used to reintroduce some of the degrees of freedom lost with LD is Dissipative Particle Dynamics (DPD). This coarse graining technique can be used to go further in the spatio-temporal scale by choosing groups of fluid particles as the simulation unit, sitting inbetween microscopic (as in MD) and macroscopic (hydrodynamic) descriptions. In practice DPD is a Langevin approach where friction acts by pairs of particles and conserves momentum.

In the standard DPD, particles interact via a soft potential, modelling the interaction between two large groups of fluid particles. In many instances DPD is used as a momentum-conserving thermostat, which thus permits to include hydrodynamics (contrary to a single Langevin approach). Local momemtum conservation results in the emergence of macroscopic hydrodynamic effects. These momentum conserving forces can then be tuned to reproduce not only thermodynamics, but also dynamical and rheological properties of diverse complex fluids. The equations of motion in DPD have the same functional form as LD and can be in fact considered as a momentum-conserving generalization of LD. The equations of motion for DPD read,

$m\vec{a} = \vec{F^c} + \vec{F^d} + \vec{F^r}.$

Where the three forces are traditionally expressed as,

$\vec{F^c}_{ij} &=\omega(r_{ij})\hat{\vec{\ppos}}_{ij}\\ \vec{F^d}_{ij} &=-\xi\omega^2(r_{ij})(\vec{\pvel}_{ij}\cdot\vec{\ppos}_{ij})\hat{\vec{\ppos}}_{ij}\\ \vec{F^r}_{ij} &=\sqrt{2\xi\kT}\omega(r_{ij})\widetilde{W}_{ij}\hat{\vec{\ppos}}_{ij}$

Where $\vec{\pvel}_{ij} = \vec{\pvel}_j - \vec{\pvel}_i$ is the relative velocity between particles $i$ and $j$ . Here $\xi$ represents a friction coefficient and is related to the random force strength via fluctuation-dissipation balance in a familiar way [2]. In general $\xi$ can be considered to be a tensorial quantity and even derived from atomistic simulations using dynamic coarse graining theory. The factor $\widetilde{W}_{ij}$ is different from the one in LD in that it affects pairs of particles (instead of each individual one), it also represents a Gaussian random number with zero mean and unit standard deviation, but must be chosen independently for each pair while ensuring symmetry so that $\widetilde{W}_{ij} = \widetilde{W}_{ji}$ . The weight function $\omega(r)$ is a soft repulsive force usually defined as

$\omega(r) = \begin{cases} \alpha\left(1-\dfrac{r}{r_{c}}\right) & r<r_{c}\\ 0 & r\ge r_{c} \end{cases}$

Where $r_{c}$ is a cut-off distance. The strength parameter, $\alpha$ , can in principle be different for each pair of particles, $i$ - $j$ , but for simplicity we will assume it is the same for every pair.

Note

The code is easily generalized for a different per-particle strength and/or friction.

Being an SDE where the forces depend on the velocities, numerical integration of the DPD equations is tricky. A simple modification can be made, sacrificing stability, by approximating the velocity to just first order in the Gronbech-Jensen update rule, so that the velocity depends only on the force for the current step. Unfortunately, this leads to artifacts in the transport properties and unacceptable temperature drifts. There are several strategies in the literature trying to overcome this, usually presented as modifications of the velocity Verlet algorithm.

All of these methods improve the accuracy of the predicted transport properties and response functions in exchange for increased computational cost. One popular approach is to simply use the energy-conserving velocity Verlet (see Molecular Dynamics) with the DPD forces. This yields “poor” stability and presents certain artifacts due to the mistreatment of the derivative of the noise term incurred by treating the DPD equations as an ordinary differential equation instead of a proper SDE. However, it is often good enough and while it might require a smaller time step to recover measurables to an acceptable tolerance it is the fastest approach and trivial to implement in a code already providing the velocity Verlet algorithm.

This is the approach used in UAMMD, where DPD is encoded as a Molecular Dynamics Integrator coupled with a Pair Forces Interactor encoding the DPD forces.

Note

The force-computing code for this module is located in the source code Interactor/Potential/DPD.cuh

Usage

A DPD Integrator is created by coupling a VerletNVE Molecular Dynamics Integrator with a DPD Potential (the Potential is supplied to a Pair Forces Interactor that can be then added to the Integrator).

The following parameters are available for the DPD Potential:

real temperature Temperature of the solvent in units of energy. This is $\kT$ in the formulas.
real cutOff The cut off, $r_c$ , for the weight function.
real gamma The friction coefficient, $\xi$ .
real A The strength of the weight function, $\alpha$ .
real dt The time step. Be sure to pass the same time step to DPD and the Integrator.

#include<uammd.cuh>
#include<Integrator/VerletNVE.cuh>
#include<Interactor/PairForces.cuh>
#include<Interactor/Potential/DPD.cuh>
using namespace uammd;
//A function that creates and returns a DPD integrator
auto createIntegratorDPD(UAMMD sim){
  Potential::DPD::Parameters par;
  par.temperature = sim.par.temperature;
  //The cut off for the weight function
  par.cutOff = sim.par.cutOff;
  //The friction coefficient
  par.gamma = sim.par.friction;
  //The strength of the weight function
  par.A = sim.par.strength;
  par.dt = sim.par.dt;
  auto dpd = make_shared<Potential::DPD>(dpd);
  //From the example in PairForces
  auto interactor = createPairForcesWithPotential(sim, dpd);
  //From the example in the MD section
  // particle velocities should not be initialized
  // by VerletNVE (initVelocities=false)
  using NVE = VerletNVE;
  NVE::Parameters params;
  params.dt = par.dt;
  params.initVelocities=false;
  verlet = make_shared<NVE>(pd,  params);
  verlet->addInteractor(interactor);
  return verlet;
}

Note

The UAMMD structure in this example is taken from the example/ folders in the repository, containing, for convenience, an instance of ParticleData and a set of parameters

Smoothed Particle Hydrodynamics

Todo

More detailed description of SPH

In SPH [3], particles are understood as interpolation points from which the fluid properties can be calculated. At the end of the day these “points” behave like particles for all intents and purposes. As a matter of fact the equations of motion for particles in SPH are just the Newton equations (see Molecular Dynamics) with a specially crafted interparticle force given by

Note

Similarly to DPD, SPH [3] is encoded as a Molecular Dynamics Integrator coupled with a special Interactor encoding the SPH forces.

The force-computing code for this module is located in the source code Interactor/SPH.cuh

$\vec{F}_i = \sum_j\left[ m_j\left(\frac{P_j}{\rho_j^2} + \frac{P_i}{\rho_i^2} + \eta_{ij}\right)\cdot\nabla_i \omega(r_{ij}/h) \right],$

where $j$ are the indices of the particles that are neighbours of $i$ (i.e. those within the support distance beyond which the interpolation kernel is zero), $m$ are the masses of the particles. Here $h$ is a length scale and $\omega(r)$ is a smooth decaying interpolation function with a close support, the so-called cubic spline kernel is commonly used,

$\omega(r) := M_4(r) = \left\{\begin{aligned} &\frac{1}{6}\left[(2-r)^3 -4(1-r)^3\right], &\quad\text{if}\quad 0\le r\le 1\\ &\frac{1}{6}(2-r)^3, &\quad\text{if}\quad 1\le r\le 2\\ &0, &\quad\text{if}\quad 0\le r> 2 \end{aligned}\right.$

Hint

The source file “Interactor/SPH/Kernel.cuh” contains the interpolation function, as well as an easy way to introduce new ones.

The density, $\rho$ , on a given particle is interpolated from its neighbours as

$\rho_i = \sum_j m_j \omega(r_{ij}/h).$

On the other hand $P$ is the pressure on the particle as given by a certain equation of state, in UAMMD’s implementation we use an ideal gas approximation given by

$P_i = K(\rho_i-\rho_0),$

Hint

The source file “Interactor/SPH.cu” contains the equation of state for the pressure, which is easily modified.

where $K$ is a gas stiffness and $\rho_0$ is a rest density.

Finally, $\eta$ is an artificial viscosity introduced to improve numerical stability, in UAMMD’s implementation we use

$\eta_{ij} = -\nu\frac{\vec{v}_{ij}\cdot \vec{r}_{ij}}{r_{ij}^2+\epsilon h^2},$

where $\nu$ is a provided viscosity, and $\vec{v}_{ij} = \vec{v}_j - \vec{v}_i$ is the relative velocity of particles i and j. $\epsilon ~ 0.001$ is introduced to prevent the singularity at $r_{ij} = 0$ .

Hint

The source file “Interactor/SPH.cu” contains the computation of the artificial viscosity, which is easily modified.

Usage

We couple an SPH Interactor with a VerletNVE Integrator.

The following parameters are available for the SPH Interactor:

real support. The length scale $h$ .
real viscosity. The prefactor for the artificial viscosity, $\mu$ .
real gasStiffness. The prefactor for the ideal gas equation of state, $K$ .
real restDensity. The rest density in the equation of state, $\rho_0$ .
Box box. A Box with the simulation domain information.

#include "Integrator/VerletNVE.cuh"
#include "Interactor/SPH.cuh"
using namespace uammd;
//SPH is handled by UAMMD as a VerletNVE integrator with a special interaction
std::shared_ptr<Integrator> createIntegratorSPH(std::shared_ptr<ParticleData> pd){
  using NVE = VerletNVE;
  NVE::Parameters par;
  par.dt = 0.1;
  par.initVelocities = false;
  auto verlet = std::make_shared<NVE>(pd, par);
  SPH::Parameters params;
  real3 L = make_real3(32,32,32);
  params.box = Box(L);
  //Pressure for a given particle "i" in SPH will be computed as gasStiffness·(density_i - restDensity)
  //Where density is computed as a function of the masses of the surroinding particles
  //Particle mass starts as 1, but you can change this in customizations.cuh
  params.support = 2.4;   //Cut off distance for the SPH kernel
  params.viscosity = 1.0;   //Environment viscosity
  params.gasStiffness = 1.0;
  params.restDensity = 1.0;
  auto sph = std::make_shared<SPH>(pd, params);
  verlet->addInteractor(sph);
  return verlet;
}

References