# Simulating Brownian Dynamics (Overdamped Langevin Dynamics) using LAMMPS

Nov 06, 2017

Update 2021.11.22. LAMMPS starts to sup­port brown­ian dy­nam­ics of­ﬁ­cially in re­cent ver­sions. See this page for the of­ﬁ­cial fix com­mand

This ar­ti­cle was orig­i­nally posted on my old Wordpress blog.

LAMMPS is a very pow­er­ful Molecular Dynamics sim­u­la­tion soft­ware I use in my daily re­search. In our re­search group, we mainly run Langevin Dynamics (LD) or Brownian Dynamics (BD) sim­u­la­tion. However, for some rea­son, LAMMPS does­n’t pro­vide a way to do Brownian Dynamics (BD) sim­u­la­tion. Both the LD and BD can be used to sam­ple cor­rect canon­i­cal en­sem­ble, which some­times also be called NVT en­sem­ble.

The BD is the large fric­tion limit of LD, where the in­er­tia is ne­glected. Thus BD is also called over­damped Langevin Dynamics. It is very im­por­tant to know the dif­fer­ence be­tween LD and BD since these two terms seems be used in­dif­fer­ently by many peo­ple which is sim­ply not cor­rect.

The equa­tion of mo­tion of LD is,

$m \ddot{\mathbf{x}} = -\nabla U(\mathbf{x}) - m\gamma \dot{\mathbf{x}}+\mathbf{R}(t)$

where $m$ is the mass of the par­ti­cle, $x$ is its po­si­tion and $\gamma$ is the damp­ing con­stant. $\mathbf{R}(t)$ is ran­dom force. The ran­dom force is sub­jected to ﬂuc­tu­a­tion-dis­si­pa­tion the­o­rem. $\langle \mathbf{R}(0)\cdot\mathbf{R}(t) \rangle = 2m\gamma\delta(t)/\beta$. $\gamma=\xi/m$ where $\xi$ is the drag co­ef­ﬁ­cient. $\mathbf{R(t)}$ is nowhere dif­fer­en­tiable, its in­te­gral is called Wiener process. Denote the wiener process as­so­ci­ated with $\mathbf{R}(t)$ as $\omega(t)$. It has the prop­erty $\omega(t+\Delta t)-\omega(t)=\sqrt{\Delta t}\theta$, $\theta$ is the Gaussian ran­dom vari­able of zero mean, vari­ance of $2m\gamma/\beta$.

$\langle \theta \rangle = 0\quad\quad\langle \theta^{2}\rangle = 2m\gamma/\beta$

The ﬁx fix langevin pro­vided in LAMMPS is the nu­mer­i­cal sim­u­la­tion of the above equa­tion. LAMMPS uses a very sim­ple in­te­gra­tion scheme. It is the Velocity-Verlet al­go­rithm where the force on a par­ti­cle in­cludes the fric­tion drag term and the noise term. Since it is just a ﬁrst or­der al­go­rithm in terms of the ran­dom noise, it can not be used for large fric­tion case. Thus the langevin fix in LAMMPS is mainly just used as a way to con­serve the tem­per­a­ture (thermostat) in the sim­u­la­tion to sam­ple the con­for­ma­tion space. However in many cases, we want to study the dy­nam­ics of our in­ter­ested sys­tem re­al­is­ti­cally where fric­tion is much larger than the in­er­tia. We need to do BD sim­u­la­tion.

For a over­damped sys­tem, $\gamma=\xi/m$ is very large, let’s take the limit $\gamma=\xi/m\to\infty$, the bath be­comes in­ﬁ­nitely dis­si­pa­tive (overdamped). Then we can ne­glect the left side of the equa­tion of LD. Thus for BD, the equa­tion of mo­tion be­comes

$\dot{\mathbf{x}}=-\frac{1}{\gamma m}\nabla U(\mathbf{x})+\frac{1}{\gamma m}\mathbf{R}(t)$

The ﬁrst or­der in­te­gra­tion scheme of the above equa­tion is called Euler-Maruyama al­go­rithm, given as

$\mathbf{x}(t+\Delta t)-\mathbf{x}(t)=-\frac{\Delta t}{m\gamma}\nabla U(\mathbf{x})+\sqrt{\frac{2\Delta t}{m\gamma\beta}}\omega(t)$

where $\omega(t)$ is the gauss­ian ran­dom vari­able with zero mean and unit vari­ance (not the same varaible ap­peared above). Since for BD, the ve­loc­i­ties are not well de­ﬁned any­more, only the po­si­tions are up­dated. The im­ple­men­ta­tion of this scheme in LAMMPS is straight­for­ward. Based on source codes fix_langevin.cpp and fix_langevin.h in the LAMMPS, I wrote a cus­tom fix of BD my­self. The core part of the code is the fol­low­ing. The whole code is on github.

void FixBD::initial_integrate(int vflag)
{
double dtfm;
double randf;

// update x of atoms in group

double **x = atom->x;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;

if (rmass) {
for (int i = 0; i < nlocal; i++)
dtfm = dtf / rmass[i];
randf = sqrt(rmass[i]) * gfactor;
x[i][0] += dtv * dtfm * (f[i][0]+randf*random->gaussian());
x[i][1] += dtv * dtfm * (f[i][1]+randf*random->gaussian());
x[i][2] += dtv * dtfm * (f[i][2]+randf*random->gaussian());
}

} else {
for (int i = 0; i < nlocal; i++)
dtfm = dtf / mass[type[i]];
randf = sqrt(mass[type[i]]) * gfactor;
x[i][0] += dtv * dtfm * (f[i][0]+randf*random->gaussian());
x[i][1] += dtv * dtfm * (f[i][1]+randf*random->gaussian());
x[i][2] += dtv * dtfm * (f[i][2]+randf*random->gaussian());
}
}
}


As one can see, the im­ple­men­ta­tion of the in­te­gra­tion scheme is easy, shown above. dtv is the time step $\Delta t$ used. dtfm is $1/(\gamma m)$ and randf is $\sqrt{2m\gamma/(\Delta t\beta)}$.

The Euler-Maruyama scheme is a sim­ple ﬁrst or­der al­go­rithm. Many stud­ies has been done on higher or­der in­te­gra­tion scheme al­low­ing large time step be­ing used. I also im­ple­mented a method shown in this pa­per. The in­te­gra­tion scheme called BAOAB is very sim­ple, given as

$\mathbf{x}(t+\Delta t)-\mathbf{x}(t)=-\frac{\Delta t}{m\gamma}\nabla U(\mathbf{x})+\sqrt{\frac{\Delta t}{2m\gamma\beta}}(\omega(t+\Delta t)+\omega(t))$

The source code of this method can be down­loaded. In ad­di­tion, feel free to fork my Github repos­i­tory for fix bd and fix bd/baoab. I have done some tests and have been us­ing this code in my re­search for a while and haven’t found prob­lems. But please test the code your­self if you in­tend to use it and wel­come any feed­back if you ﬁnd any prob­lems.

To de­cide whether to use LD or BD in the sim­u­la­tion, one need to com­pare rel­e­vant timescales. Consider a free par­ti­cle gov­erned by the Langevin equa­tion. Solving for the ve­loc­ity au­to­cor­re­la­tion func­tion leads to, $\langle v(0)v(t)\ran­gle=(kT/​m)e^{-\gamma t}$. This shows that the re­lax­ation time for mo­men­tum is $\tau_m = 1/\gamma=m/\xi$. There is an­other timescale called Brownian timescale cal­cu­lated by $\tau_{BD}=\sigma^2\xi/kT$ where $\sigma$ is the size of the par­ti­cle. $\tau_{BD}$ is the timescale at which the par­ti­cle dif­fuses about its own size. If $\tau_{BD}\gg \tau_m$ and if you are not in­ter­ested at the dy­nam­ics on the timescale $\tau_m$, then one can use BD since the mo­men­tum can be safely in­te­grated out. However, if these two timescales are com­pa­ra­ble or $\tau_{BD}<\tau_m$, then only LD can be used be­cause the mo­men­tum can­not be ne­glected in this case. To make the prob­lem more com­pli­cated, there are more than just these two timescales in most of sim­u­la­tion cases, such as the re­lax­ation time of bond vi­bra­tion, etc… Fortunately, prac­ti­cally, com­par­ing these two timescales is good enough for many cases.