Implementation of novel boundary conditions in LAMMPS

Dear all,

I am performing dissipative particle dynamics (DPD) simulations of systems comprising polymer, solvent and solid walls of arbitrarily complex geometries. The particles in my system are classified into three types: polymer, solvent and solid. The particles in the last type comprising the solid mass are frozen in place.

I wish to implement the solid-fluid boundary conditions (BCs) described in pages 118 and 119 of the following reference in LAMMPS:

Zhang, D., Shangguan, Q. & Wang, Y. An easy-to-use boundary condition in dissipative particle dynamics system. Comput. Fluids 166, 117–122 (2018).

The modified velocity verlet (vv) algorithm described therein involves the following four steps:

  1. Update the positions of all particles
  2. Determine the predicted velocity of each particle
  3. Calculate forces based on the predicted velocity
  4. Update the velocities of all particles

To implement these BCs, I need to determine the minimum distance, r_min, between each mobile (polymer or solvent) bead and the solid beads within the cutoff radius,r_c, at the end of step 1. I also need to determine the total force, F_s, exerted on each mobile bead by the solid beads within r_c. Finally, I need to modify the velocities of mobile beads determined in step 4 based on the values of F_s and r_min. Could someone kindly advise as to the most efficient means of accomplishing this using the 11Aug17 version of LAMMPS? In particular, can it be done using existing functionality or will an additional fix be required?

Best,
Karthik.

***DISCLAIMER*** The sender of this email is an alumnus of National University of Singapore (NUS). Kindly note that NUS is not responsible for the contents of this email, and views and opinions expressed are solely the sender's.

I think you would need to write a new time-integration fix. You can

use one of the many existing ones as a starting template.

Steve

In my case, velocity verlet integration proceeds as per normal. The only difference is that there is an additional step which involves correcting the velocities of those particles for which the velocity, r_min and F_s satisfy certain conditions. Therefore, would it be sufficient to write a fix that performs this correction at the end of every time-step? Also, can I determine r_min and F_s in the same loop that performs pairwise force calculations in pair_dpd.cpp?

Best,

Karthik.

can I determine r_min and F_s in the same loop that performs pairwise force calculations in pair_dpd.cpp?

Not unless you want to edit pair dpd, which means we likely won’t accept your changes

for inclusion in LAMMPS. But maybe you don’t care about that.

Fixes can requests neighbor lists (or copies of an existing one, so no overhead to build them).

They can loop over the list and calculate things. I suggest you just program up a new

fix that does what you want and not worry about its overlap with other parts of LAMMPS.

Steve

I tried to follow your suggestion of putting all the relevant code in a single fix without editing pair_dpd. I called the fix test/bc1. The contents of ‘fix_test_bc1.h’ and ‘fix_test_bc1.cpp’ are as follows:

fix_test_bc1.h

#ifdef FIX_CLASS

FixStyle(test/bc1,FixTestBc1)

#else

#ifndef LMP_FIX_TEST_BC1_H
#define LMP_FIX_TEST_BC1_H

#include “fix.h”

namespace LAMMPS_NS {

class FixTestBc1 : public Fix {

friend class PairDPD;

public:
FixTestBc1(class LAMMPS *, int, char **);
virtual ~FixTestBc1();
int setmask();
void init();
virtual void post_force();
virtual void end_of_step();
double memory_usage();

protected:

double alphaone,alphatwo ;
double *rminsolid,**fsolid;
int maxatom;

};

}

#endif
#endif

fix_test_bc1.cpp

#include <math.h>
#include <string.h>
#include <stdlib.h>
#include “fix_test_bc1.h”
#include “math_const.h”
#include “atom.h”
#include “atom_vec.h”
#include “comm.h”
#include “update.h”
#include “force.h”
#include “input.h”
#include “neighbor.h”
#include “neigh_list.h”
#include “neigh_request.h”
#include “random_mars.h”
#include “memory.h”
#include “error.h”
#include “pair_dpd.h”

using namespace LAMMPS_NS;

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

FixTestBc1::FixTestBc1(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)

{
if (narg < 6) error->all(FLERR,“Illegal fix test/bc1 command”);

nevery = force->inumeric(FLERR,arg[3]);
if (nevery != 1) error->all(FLERR,“Illegal fix test/bc1 command”);
alphaone = force->numeric(FLERR,arg[4]);
alphatwo = force->numeric(FLERR,arg[5]);

maxatom = 1;

memory->create(rminsolid,maxatom,“test/bc1:rminsolid”);
memory->create(fsolid,maxatom,3,“test/bc1:fsolid”);

}

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

FixTestBc1::~FixTestBc1()
{
memory->destroy(rminsolid);
memory->destroy(fsolid);
}

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

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

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

void FixTestBc1::init()

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

}

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

void FixTestBc1::post_force()

{

int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,randnum,factor_dpd;
int *ilist,*jlist,*numneigh,**firstneigh;

double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int n = atom->ntypes; //additions
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
double dtinvsqrt = 1.0/sqrt(update->dt);

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

maxatom = atom->nmax;
memory->destroy(fsolid);
memory->destroy(rminsolid);
memory->create(rminsolid,maxatom,“test/bc1:rminsolid”);
memory->create(fsolid,maxatom,3,“test/bc1:fsolid”);

for (ii = 0; ii < maxatom; ii++) {
rminsolid[ii]=1000.0;
for (jj=0; jj <= 2; jj++) {
fsolid[ii][jj] = 0.0;
}
}

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itype = type[i];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];

jlist = firstneigh[i];
jnum = numneigh[i];

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

if (rsq < cutsq[itype][jtype]) {
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 = delxdelvx + delydelvy + delz*delvz;
wd = 1.0 - r/cut[itype][jtype];
randnum = random->gaussian();

// conservative force = a0 * wd
// drag force = -gamma * wd^2 * (delx dot delv) / r
// random force = sigma * wd * rnd * dtinvsqrt;

fpair = a0[itype][jtype]wd;
fpair -= gamma[itype][jtype]wdwd
dotrinv;
fpair += sigma[itype][jtype]wdrandnum
dtinvsqrt;
fpair = factor_dpdrinv;

if (jtype==n) {

fsolid[i][0] += delxfpair;
fsolid[i][1] += dely
fpair;
fsolid[i][2] += delz*fpair;
if (r < rminsolid[i]) rminsolid[i] = r ;

}

if (newton_pair || j < nlocal) {
if (itype==n){
fsolid[j][0] -= delxfpair;
fsolid[j][1] -= dely
fpair;
fsolid[j][2] -= delz*fpair;
if (r < rminsolid[j]) rminsolid[j] = r ;
}
}

}
}
}

}

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

void FixTestBc1::end_of_step()

{
int i,itype;
int n = atom->ntypes;
int *type = atom->type;
int nlocal = atom->nlocal;
double **v = atom->v;
double **f = atom->f;

double r,fone,dot,fmag,fnorm[3],vmag,beta,ftwo ;

for (i = 0; i < nlocal; i++) {
itype = type[i];
if (itype!=ntype) {
r=rminsolid[i];
if (r < cut_global) {
fone = pow((cut_global-r)/cut_global,alphaone);
dot = fsolid[i][0]v[i][0] + fsolid[i][1]v[i][1] + fsolid[i][2]v[i][2];
if (dot <= 0) {
fmag=sqrt(fsolid[i][0]fsolid[i][0] + fsolid[i][1]fsolid[i][1] + fsolid[i][2]fsolid[i][2]);
fnorm[0]=fsolid[i][0]/fmag;
fnorm[1]=fsolid[i][1]/fmag;
fnorm[2]=fsolid[i][2]/fmag;
vmag = sqrt( v[i][0]v[i][0] + v[i][1]v[i][1] + v[i][2]v[i][2]);
dot /= (-1.0 * vmag * fmag);
beta=acos(dot) * 180.0 / MY_PI;
ftwo = pow((90.0-beta)/90.0,alphatwo);
v[i][0] += vmag
fone
ftwo
fnorm[0];
v[i][1] += vmag
fone
ftwo
fnorm[1];
v[i][2] += vmag
fone
ftwo
fnorm[2];
}
}

}
}
}

/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */

double FixTestBc1::memory_usage()
{

double bytes=0.0;
bytes += maxatomsizeof(double);
bytes += maxatom
3*sizeof(double);

return bytes;
}

When I attempt to compile LAMMPS, I get the following error messages:

…/fix_test_bc1.cpp: In member function ‘virtual int LAMMPS_NS::FixTestBc1::setmask()’:
…/fix_test_bc1.cpp:73: error: ‘POST_FORCE’ was not declared in this scope
…/fix_test_bc1.cpp:74: error: ‘END_OF_STEP’ was not declared in this scope
…/fix_test_bc1.cpp: In member function ‘virtual void LAMMPS_NS::FixTestBc1::post_force()’:
…/fix_test_bc1.cpp:116: error: ‘list’ was not declared in this scope
…/fix_test_bc1.cpp:151: error: ‘sbmask’ was not declared in this scope
…/fix_test_bc1.cpp:160: error: ‘cutsq’ was not declared in this scope
…/fix_test_bc1.cpp:162: error: ‘EPSILON’ was not declared in this scope
…/fix_test_bc1.cpp:168: error: ‘cut’ was not declared in this scope
…/fix_test_bc1.cpp:169: error: request for member ‘gaussian’ in ‘random’, which is of non-class type ‘long int()throw ()’
…/fix_test_bc1.cpp:175: error: ‘a0’ was not declared in this scope
…/fix_test_bc1.cpp:176: warning: pointer to a function used in arithmetic
…/fix_test_bc1.cpp:176: warning: pointer to a function used in arithmetic
…/fix_test_bc1.cpp:176: error: invalid operands of types ‘double(double)’ and ‘double’ to binary ‘operator*’
…/fix_test_bc1.cpp:177: error: ‘sigma’ was not declared in this scope
…/fix_test_bc1.cpp: In member function ‘virtual void LAMMPS_NS::FixTestBc1::end_of_step()’:
…/fix_test_bc1.cpp:226: error: ‘ntype’ was not declared in this scope
…/fix_test_bc1.cpp:228: error: ‘cut_global’ was not declared in this scope
…/fix_test_bc1.cpp:238: error: ‘MY_PI’ was not declared in this scope

May I know where the issue lies? I need to access some variables such as the cutoff radius for pairwise interactions and the DPD force coefficients defined and used in the pair_dpd files but I can’t seem to access them.

DISCLAIMER The sender of this email is an alumnus of National University of Singapore (NUS). Kindly note that NUS is not responsible for the contents of this email, and views and opinions expressed are solely the sender’s.

I tried to follow your suggestion of putting all the relevant code in a
single fix without editing pair_dpd. I called the fix test/bc1. The
contents of 'fix_test_bc1.h' and 'fix_test_bc1.cpp' are as follows:

​[...]


#include "pair_dpd.h"

using namespace LAMMPS_NS;

​there should be a:

using namespace FixConst;

here:​

[
​...]

   mask |= POST_FORCE;

   mask |= END_OF_STEP;

or ​you need to use

FixConst::POST_FORCE

and

FixConst::END_OF_STEP

here.​

​axel.

Dear Axel and Steve,

Thank you for your suggestions. I am facing a slight issue. What I want to do is to modify the velocities of the mobile atoms in my system at the end of each time-step. However, in order to do this, I need to know the force exerted upon each mobile bead by all the solid (frozen) beads within the cutoff radius, which is calculated earlier in the time-step using the code in pair_dpd. In DPD, the pairwise forces have a random component. Therefore, if I write my own fix that repeats the force calculation in pair_dpd, the random force that I get will be different from that calculated by pair_dpd earlier in the time-step. The resulting dynamics will be different from the one obtained using the method of the reference in my earlier message.

Hence, it seems that I have no choice but to edit pair_dpd. If I define additional arrays in pair_dpd, how can I ensure that their contents will be available for my fix to use as input at the end of each time-step?

DISCLAIMER The sender of this email is an alumnus of National University of Singapore (NUS). Kindly note that NUS is not responsible for the contents of this email, and views and opinions expressed are solely the sender’s.

Pair styles support an extract() method that other classes (e.g. a Fix) can use

to request a pointer to info stored by the pair style. Grep existing pair

styles for extract and you will see some examples.

Steve