Simulating Brownian Dynamics (Overdamped Langevin Dynamics) using LAMMPS

Update 2021.11.22. LAMMPS starts to sup­port brown­ian dy­nam­ics of­fi­cially in re­cent ver­sions. See this page for the of­fi­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,

mx¨=U(x)mγx˙+R(t)m \ddot{\mathbf{x}} = -\nabla U(\mathbf{x}) - m\gamma \dot{\mathbf{x}}+\mathbf{R}(t)

where mm is the mass of the par­ti­cle, xx is its po­si­tion and γ\gamma is the damp­ing con­stant. R(t)\mathbf{R}(t) is ran­dom force. The ran­dom force is sub­jected to fluc­tu­a­tion-dis­si­pa­tion the­o­rem. R(0)R(t)=2mγδ(t)/β\langle \mathbf{R}(0)\cdot\mathbf{R}(t) \rangle = 2m\gamma\delta(t)/\beta. γ=ξ/m\gamma=\xi/m where ξ\xi is the drag co­ef­fi­cient. R(t)\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 R(t)\mathbf{R}(t) as ω(t)\omega(t). It has the prop­erty ω(t+Δt)ω(t)=Δtθ\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γ/β2m\gamma/\beta.

θ=0θ2=2mγ/β\langle \theta \rangle = 0\quad\quad\langle \theta^{2}\rangle = 2m\gamma/\beta

The fix 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 first 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, γ=ξ/m\gamma=\xi/m is very large, let’s take the limit γ=ξ/m\gamma=\xi/m\to\infty, the bath be­comes in­fi­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

x˙=1γmU(x)+1γmR(t)\dot{\mathbf{x}}=-\frac{1}{\gamma m}\nabla U(\mathbf{x})+\frac{1}{\gamma m}\mathbf{R}(t)

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

x(t+Δt)x(t)=ΔtmγU(x)+2Δtmγβω(t)\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 ω(t)\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­fined 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 *mask = atom->mask;
  int nlocal = atom->nlocal;
  if (igroup == atom->firstgroup) nlocal = atom->nfirst;

  if (rmass) {
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit) {
        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++)
      if (mask[i] & groupbit) {
        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 Δt\Delta t used. dtfm is 1/(γm)1/(\gamma m) and randf is 2mγ/(Δtβ)\sqrt{2m\gamma/(\Delta t\beta)}.

The Euler-Maruyama scheme is a sim­ple first 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

x(t+Δt)x(t)=ΔtmγU(x)+Δt2mγβ(ω(t+Δt)+ω(t))\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 find 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, v(0)v(t)=(kT/m)eγt\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 τm=1/γ=m/ξ\tau_m = 1/\gamma=m/\xi. There is an­other timescale called Brownian timescale cal­cu­lated by τBD=σ2ξ/kT\tau_{BD}=\sigma^2\xi/kT where σ\sigma is the size of the par­ti­cle. τBD\tau_{BD} is the timescale at which the par­ti­cle dif­fuses about its own size. If τBDτm\tau_{BD}\gg \tau_m and if you are not in­ter­ested at the dy­nam­ics on the timescale τm\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 τBD<τm\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.