Simulating Brownian Dynamics (Overdamped Langevin Dynamics) using LAMMPS

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 fluc­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­fi­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 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, $\gamma=\xi/m$ is very large, let’s take the limit $\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

$\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

$\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­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 $\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 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

$\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.