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.7. How a timestep works — LAMMPS documentation

Got it. Thank you. :pray:

Dear Dr. Kohlmeyer,

I followed your suggestion and wrote the fix, however the momentum is still not conserved. Please see the attached code for the fix. I understand and don’t expect you to debug my code, but I request you to see if I implemented your suggestions correctly.

#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 <iostream>
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), vmod(nullptr), randla(nullptr)
{
	if (narg != 7) error->all(FLERR,"Illegal fix temp/lowe command");

	// Lowe-Andersen thermostat should be applied every step
	nevery = 1;
    global_freq = nevery;

	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

	randla = new RanMars(lmp, seedla);

    peratom_flag = 1;
    size_peratom_cols = 3;
    peratom_freq = 1;

    nmax = 0; // Initialization 
 
	comm_reverse = 3;
}
/* ---------------------------------------------------------------------- */
FixTempLowe::~FixTempLowe()
{
	delete randla;
	memory->destroy(vmod);

}
/* ---------------------------------------------------------------------- */
void FixTempLowe::init()
{
    if (force->newton_pair == 0) error->all(FLERR, "Fix temp/lowe requires newton pair on");
	// needs half neighbor list
	neighbor->add_request(this);

}
/* ---------------------------------------------------------------------- */
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()
{
    if ( (atom->nmax) > nmax ) {
        memory->destroy(vmod);
        nmax = atom->nmax;
        memory->create(vmod, nmax, 3, "temp/lowe:vmod");
        array_atom = vmod;
    
        // Initialize
        for (int n = 0; n<nmax; n++) {
            vmod[n][0] = vmod[n][1] = vmod[n][2] = 0.0;
	    }
    }

	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,delvij;
	double rsq,r,rinv,lambda,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;
    int nghost = atom->nghost;
    int nall = nlocal+nghost;
    int newton_pair = force->newton_pair;

    // Copy of atom velocities
    for (int k = 0; k<nall; k++) {
        vmod[k][0] = v[k][0];
        vmod[k][1] = v[k][1];
        vmod[k][2] = v[k][2];
    }
 
	inum = list->inum;
	ilist = list->ilist;
	numneigh = list->numneigh;
	firstneigh = list->firstneigh;
    
	// loop over neighbors of my atoms
    double rcutsq = rcutla*rcutla;
	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 = vmod[i][0];
		vytmp = vmod[i][1];
		vztmp = vmod[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 >= (randla->uniform())) ) {
				r = sqrt(rsq);
				if (r < EPSILON) continue; // r can be 0.0 in DPD systems
				rinv = 1.0/r;
				delvx = vxtmp - vmod[j][0];
				delvy = vytmp - vmod[j][1];
				delvz = vztmp - vmod[j][2];
				dotij = (delx*delvx + dely*delvy + delz*delvz)*rinv;

				muij = mass[itype]*mass[jtype]/(mass[itype]+mass[jtype]);
				sigmaij = sqrt((force->boltz)*templa/muij);
				lambda = sigmaij*(randla->gaussian());
				delvij = (lambda-dotij)*rinv;

                // Relative velocity
				vmod[i][0] += delvij*delx*muij/mass[itype];
				vmod[i][1] += delvij*dely*muij/mass[itype];
				vmod[i][2] += delvij*delz*muij/mass[itype];
                
                vmod[j][0] -= delvij*delx*muij/mass[jtype];
				vmod[j][1] -= delvij*dely*muij/mass[jtype];
				vmod[j][2] -= delvij*delz*muij/mass[jtype];
                
			} // If condition 
  
		} // For over jj

	} // For over ii
  
    // Assign modified velocities to local atoms
    for (int k = 0; k < nlocal; k++) {
        v[k][0] = vmod[k][0]; 
        v[k][1] = vmod[k][1]; 
        v[k][2] = vmod[k][2];
    }

    // Communicate modified velocities to ghost atoms
	if (newton_pair) comm->reverse_comm(this);
    
} // End
/* ---------------------------------------------------------------------- */
int FixTempLowe::pack_reverse_comm(int n, int first, double *buf)
{
	int i, m, last;

	m = 0;
  	last = first + n;
    for (i = first; i < last; i++) {
        buf[m++] = vmod[i][0];
      	buf[m++] = vmod[i][1];
      	buf[m++] = vmod[i][2];
    }
  	return m;
}
/* ---------------------------------------------------------------------- */
void FixTempLowe::unpack_reverse_comm(int n, int *list, double *buf)
{
  	int i, j, m;
    double **v = atom->v;
  	
    m = 0;
  	for (i = 0; i < n; i++) {
      	j = list[i]; 
      	v[j][0] = buf[m++];
      	v[j][1] = buf[m++];
      	v[j][2] = buf[m++];
    }
}
/* ---------------------------------------------------------------------- */
double FixTempLowe::memory_usage()
{
    double bytes = (double) (atom->nmax)*3*sizeof(double);
    return bytes;
}

Regards
VSA

How are those two things different?

Have you also created a custom DPD pair style that does not include its own thermostat? Otherwise you would be using a pairwise thermostat twice!

Have you checked with the no-thermostat-DPD model that your momentum is conserved with fix nve?

Have you checked your thermostat on a simple example, e.g. the LJ melt example to check of you get proper thermostatting and the expected momentum conservation there. For example by comparing with fix nve and fix nve + fix langevin or fix nvt. Of course, you need to extend the duration of the run a bit and separate out the equilibration. This would be standard practice for testing/debugging a new thermostat style.

I am not creating a custom pair style, rather I simply set the \gamma = 0 using the pair_coeff. Doing so will not maintain the temperature at a desired value but momentum is conserved.

Actually, I am confident of troubleshooting the physics part of it, what I am not fully confident of is the coding part (packing and unpacking functions, flags etc). I tried my best to go thru multiple fixes to code but just to be sure I wanted an expert’s assistance. Further, I successfully implemented the thermostat in my serial code.

Regards
VSA

But that can only be ascertained by creating some simple test inputs and then setting suitable breakpoints and then following the code in a debugger (or do non-debugger approach of peppering the code with print statements). That is debugging. If anybody takes a (superficial) look at the code and says “it looks correct”, that does not confirm its correctness. You can emulate the necessary communication pattern with a single MPI ranks, where pairs of atoms are interacting across periodic boundaries. Then one of the two will be a ghost atom and its added velocity would have to be added to its original (a local atom with the same atom ID) with a reverse communication.

Thank you Dr. Kohlmeyer for your suggestions. It’s kind of working, the momentum is conserved now. However, I am a bit uncomfortable with the way I got it to work. Here are the steps I followed on a single MPI ranks.

  1. Create a 2d-box of 4x4 size with 2 atoms (one at (-3.8,-3.8) and the other at (3.8,3.8)) with periodic boundary conditions and rcut = 1.

  2. Copy the original velocities (A total of eight entries, 2-local, 6-ghost) to a local variable, vmod.

  3. Modify vmod variable as follows
    3.1. Store the modified velocities in the vmod as I loop over the pairs. In the present scenario, there is only one pair. Therefore only two entries in the vmod are modified.
    3.2 Before performing a reverse communication, go over and correct the unmodified entries in vmod based on global ID. (This is the step I am uncomfortable with)

  4. Perform a reverse communication with pack and unpack functions.

Regards
VSA

From your description and a quick cross-check with your code, I can say with confidence that you have neither followed my suggestions, nor properly understood what reverse and forward communications do and what they are used for, nor actually applied them correctly. You should see in your debugging that your code is not doing what your algorithm states it should.

I have given this all the attention I want to give, provided all the information or pointers to it that would be necessary to implement what you need. As I now stated multiple times: the remaining work is just careful re-reading my suggestions, plain debugging and logical thinking. I have no intention to do that for you.

Seems to me like this Lowe-Andersen thermostat is basically a DPD thermostat with the following changes:

  • velocity thermostatting activates stochastically instead of for every pair in cutoff distance
  • no ramp function (or in the LAMMPS documentation notation – w(r) is 1 within the cutoff and 0 outside)
  • \gamma is calculated per-pair to completely cancel the previous relative velocity instead of partially cancelling it
  • \sigma is rescaled appropriately (the documentation formula might even be correct with the redefined \gamma).

So why not just copy the pair dpd or pair dpd/tstat styles and make those changes? You can implement the velocity update as a force change and use the fix mvv/tstat integrator for fine-tuning the equations of motion (it seems like \lambda = 1 in the integrator should get you what you need). Then you can use all of LAMMPS’s existing comm infrastructure for free – including the per-atom virial calculations.

Thank you Dr. Kohlmeyer and srtee. Now, I got it working.