Lowe-Andersen Thermostat in LAMMPS

Dear LAMMPS community,

I am trying to implement Lowe-Andersen thermostat for a DPD system. I have written a fix that’s kind of working but not quite right. I believe the issue is with communicating the modified velocity between pairs.


Brief overview of Lowe-Andersen thermostat

After the Verlet integration is completed, a pair of atoms (within the cutoff) are selected with a probability and their relative velocity is exchanged with a reservoir following M-B statistics. The probability = freq*dt, where freq is set by the user and dt is the timestep.

I am not sure if this is the best way, but I request a full neighbor list so that inter-processor communication is avoided. The issue is the temperature of the system constantly drifts to the higher side. This, I believe is due to velocities of both the atoms (in a given pair) are not getting modified by exactly equal and opposite amount (hence conservation of momentum). I am thinking it might be an issue with the ghost atoms.

LAMMPS version: parallel build, 29Sep2021

=========
lammps inscript to reproduce the observed behavior:

dimension 2
units lj
comm_modify vel yes
newton on

atom_style atomic
atom_modify map yes
neighbor 0.3 bin
neigh_modify delay 0 every 1 check yes

region domain block -5.0 5.0 -5.0 5.0 -0.1 0.1 units box
create_box 1 domain

variable seed equal ‘round(random(1, 899999999, 868404421))’

create_atoms 1 random 300 $(v_seed) domain

variable usrtemp equal 1.0
variable usrrc equal 1.0
variable usrfreq equal 10.0

mass * 1.0

pair_style dpd {usrtemp} {usrrc} $(v_seed)

pair_coeff * * 25.0 0.0

thermo 1000
timestep 0.01

velocity all create {usrtemp} (v_seed) dist gaussian

fix f4 all nve
fix f3 all temp/lowe {usrtemp} {usrrc} {usrfreq} (v_seed)
fix f5 all enforce2d

run 10000

=============
fix_temp_lowe.h

#ifdef FIX_CLASS
// clang-format off
FixStyle(temp/lowe,FixTempLowe);
// clang-format on
#else

#ifndef LMP_FIX_TEMP_LOWE_H
#define LMP_FIX_TEMP_LOWE_H

#include “fix.h”

namespace LAMMPS_NS {

class FixTempLowe : public Fix {
public:
FixTempLowe(class LAMMPS *, int, char **);
int setmask();
void init();
void init_list(int, class NeighList *);
void end_of_step();

private:
double freqLA;
double rcutLA;
double tempLA;
int seedLA;
class RanMars *random;
int me;

class NeighList *list;

};

} // namespace LAMMPS_NS

#endif
#endif

===============
fix_temp_lowe.cpp

#include “fix_temp_lowe.h”

#include “atom.h”
#include “error.h”
#include “force.h”
#include “memory.h”
#include “neigh_list.h”
#include “neigh_request.h”
#include “neighbor.h”
#include “pair.h”
#include “update.h”

#include “random_mars.h”
#include “comm.h”

#include
#include

#include
using namespace std;

using namespace LAMMPS_NS;
using namespace FixConst;

#define EPSILON 1.0e-10

/* ---------------------------------------------------------------------- */

FixTempLowe::FixTempLowe(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), list(nullptr)
{
if (narg != 7) error->all(FLERR,“Illegal fix temp/lowe command”);

// Lowe-Andersen thermostat should be applied every step

nevery = 1;

tempLA = utils::numeric(FLERR,arg[3],false,lmp); // set temperature
rcutLA = utils::numeric(FLERR,arg[4],false,lmp); // set cut-off
freqLA = utils::numeric(FLERR,arg[5],false,lmp); // set frequency
seedLA = utils::inumeric(FLERR,arg[6],false,lmp); // seed

random = new RanMars(lmp, seedLA + comm->me);

}
/* ---------------------------------------------------------------------- */

void FixTempLowe::init()
{
// needs full neighbor list

int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}

/* ---------------------------------------------------------------------- */

void FixTempLowe::init_list(int /id/, NeighList ptr)
{
list = ptr;
}
/
---------------------------------------------------------------------- */

int FixTempLowe::setmask()
{
int mask = 0;
mask |= END_OF_STEP;
return mask;
}

/* ---------------------------------------------------------------------- */

void FixTempLowe::end_of_step()
{

int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double muij,sigmaij,delvi,delvj;
double rsq,r,rinv,dot,dotij;
int *ilist,*jlist,*numneigh,**firstneigh;

double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;

inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;

double rcutsq = rcutLArcutLA;
double prob = freqLA
(update->dt);

// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
  j = jlist[jj];
  j &= NEIGHMASK;
  delx = xtmp - x[j][0];
  dely = ytmp - x[j][1];
  delz = ztmp - x[j][2];
  rsq = delx*delx + dely*dely + delz*delz;
  jtype = type[j];

  if ( (rsq < rcutsq) && (prob >= (random->uniform())) ) {
r = sqrt(rsq);
if (r < EPSILON) continue; // r can be 0.0 in DPD systems
rinv = 1.0/r;
    delvx = vxtmp - v[j][0];
    delvy = vytmp - v[j][1];
    delvz = vztmp - v[j][2];
    dot = (delx*delvx + dely*delvy + delz*delvz)*rinv;
   
muij = mass[type[i]]*mass[type[j]]/(mass[type[i]]+mass[type[j]]);
sigmaij = sqrt((force->boltz)*tempLA/muij);
dotij = sigmaij*(random->gaussian());
delvi = muij/mass[type[i]]*(dotij-dot)*rinv;
delvj = muij/mass[type[j]]*(dotij-dot)*rinv;

v[i][0] += delvi*delx;
    v[i][1] += delvi*dely;
    v[i][2] += delvi*delz;

// Not sure about the following
// if (j < nlocal) {
v[j][0] -= delvjdelx;
v[j][1] -= delvj
dely;
v[j][2] -= delvj*delz;
// }

  } // If condition 
  
} // For over jj

} // For over ii

} // End

Thanks and regards
Sathish

Why not use pair style dpd/tstat?

Dear Dr. Kohlmeyer,

Thank you for your response.

The DPD has an in-built thermostat, but the corresponding hydrodynamics are gas-like (very low Schmidt number). The Lowe-Andersen thermostat yields better hydrodynamics (Schmidt number ~ 1000). By tuning the frequency of velocity exchange with the reservoir, one has flexibility over Schmidt number. These properties are absent in a dpd-type thermostat.

Thanks and regards
Sathish

If you are developing new code, you should be doing this with the latest development version on a feature branch via git. Otherwise, you run a high risk that your code will be obsolete before it is done, the online developer documentation not fully applicable, and extra porting effort required to integrate it into the LAMMPS distribution.

If you want the velocity changes per pair of atoms to be consistent, then you cannot avoid communication. Here is how you can do this with the internal features of LAMMPS:

  • request a half neighbor list and ensure that the newton pair setting is on
  • signal that you intend to do a reverse communication for 3 items
  • add pack/unpack functions for the reverse communication that packs/unpacks the “vmod” data (see below)
  • set up an array double **vmod with the dimensions atom->nmax and 3, delete and reallocate if atom->nmax grows.
  • set this array to zero
  • process the neighbor list and store velocity changes in vmod for both atoms of a pair
  • do a reverse communication to add the content of vmod from ghost atoms to their respective local atoms
  • apply the combined velocity changes to the atom velocities (for local atoms only)

Thank you Dr. Kohlmeyer.

Is there any fix where this kind of inter processor communication happens? I will start from there.

Please first read 4.4.2. Communication — LAMMPS documentation and
4.5. Communication patterns — LAMMPS documentation

Fix styles that use these functions can be easily found using the grep command.

I am working on your suggestions. Thank you.

Excuse me for my naive question. Why this communication is not an issue with forces? For example (from paid_dpd.cpp), the following lines are modifying the forces of pairs of particles.

    fpair = a0[itype][jtype]*wd;
    fpair -= gamma[itype][jtype]*wd*wd*dot*rinv;
    fpair += sigma[itype][jtype]*wd*randnum*dtinvsqrt;
    fpair *= factor_dpd*rinv;

    f[i][0] += delx*fpair;
    f[i][1] += dely*fpair;
    f[i][2] += delz*fpair;
    if (newton_pair || j < nlocal) {
      f[j][0] -= delx*fpair;
      f[j][1] -= dely*fpair;
      f[j][2] -= delz*fpair;
    }

It is. But since this applies to all force computations, the corresponding reverse communication is done as part of the time integration class and thus a lot of redundant code avoided. Please see the workflow pseudo code on this page of the documentation: 4.6. How a timestep works — LAMMPS documentation

Got it. Thank you. :pray: