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] -= delvjdely;
v[j][2] -= delvj*delz;
// }
} // If condition
} // For over jj
} // For over ii
} // End
Thanks and regards
Sathish
