Trying to change fix viscous so that it acts perpendicularly to a graphene plane

Hi all,

I am trying to implement a change in fix viscous so that the force acts only perpendicularly to my graphene sheet. The code can be found at this link in the first appendix:

I am using cmake to recompile lammps, and can successfully compile it with the default code found in fix viscous:

void FixViscous::post_force(int /*vflag*/)
  // apply drag force to atoms in group
  // direction is opposed to velocity vector
  // magnitude depends on atom type

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

  double drag;

  for (int i = 0; i < nlocal; i++)
    if (mask[i] & groupbit) {
      drag = gamma[type[i]];
      f[i][0] -= drag*v[i][0];
      f[i][1] -= drag*v[i][1];
      f[i][2] -= drag*v[i][2];

BUT, when I change it to this code in the paper above:

void FixViscous::post_force(int /*vflag*/)
  double rssq = 25; //square of rs
  double xtmp, ytmp, ztmp, delx, dely, delz, rsq;
  double xreg[100], yreg[100], zreg[100];
  double xx, yy, zz, xy, xz, yz, xsum, ysum, zsum, vdotn;
  double A[3][3], b[3], n[3], nmag;
  int inum, jnum, i, ii, j, jj, nreg, k;
  int *ilist,*jlist,*numneigh,**firstneigh;
  double shift = 1.03040961; //arbitrary shift to prevent illconditioning
  double **x = atom->x;
  double **v = atom->v;
  double **f = atom->f;
  int *mask = atom->mask;
  int *type = atom->type;
  int nlocal = atom->nlocal;
  // invoke full neighbor list (will copy or build if necessary)
  // neighbor->build_one(list);
  inum = list->inum;
  ilist = list->ilist;
  numneigh = list->numneigh;
  firstneigh = list->firstneigh;
  double drag;
  // determine the local normal direction
  // based on all neighbors within rs of
  // the atom, then impose the drag force in
  // that direction
  for (ii = 0; ii < inum; ii++) {
    i = ilist[ii];
    nreg = 0;
    if (mask[i] & groupbit) {
      xtmp = x[i][0];
      ytmp = x[i][1];
      ztmp = x[i][2];
      jlist = firstneigh[i];
      jnum = numneigh[i];
      //add point to regression array
      xreg[nreg] = xtmp + shift;
      yreg[nreg] = ytmp + shift;
      zreg[nreg] = ztmp + shift;
      nreg += 1;
      for (jj = 0; jj < jnum; jj++) {
        j = ilist[jj];
        if (mask[j] & groupbit) {
          delx = xtmp - x[j][0];
          dely = ytmp - x[j][1];
          delz = ztmp - x[j][2];
          domain -> minimum_image(delx,dely,delz);
          rsq = delx*delx + dely*dely + delz*delz;
          if (rsq < rssq) {
            xreg[nreg] = x[j][0] + shift;
            yreg[nreg] = x[j][1] + shift;
            zreg[nreg] = x[j][2] + shift;
            nreg += 1;

the compiling fails and I get this error:

What is going on here?

If I uncomment this line: // neighbor->build_one(list);
then the error becomes:

Do you what is going on here/how to fix this??

This appears to be very incomplete.

Just standard C++ stuff, you are using a variable (list) that has not been properly declared.
And you are using a class (Domain) where you only have a forward declaration since you apparently didn’t add the corresponding include file.

This is also a standard C++ thing, where you are using a member function of a class where you have only a forward declaration and are missing the include file declaring the class.
Also uncommenting the build_one() call is incorrect. You would need to request a perpetual, full neighbor list elsewhere in the class implementation and then the neighbor list code would figure out the most effective way to provide it while it generates also the neighbor list for the pair style.

Overall, there is a lot missing that can only be guessed. As a minimum, there needs to be a “list” pointer declared in the class header file as a (protected) class member. You need to add a neighbor list request and the fix needs to have a list_init() function so that the (perpetual) neighbor list will be assigned to the “list” pointer after the neighbor list code has built it.

There are likely more things missing, but that is very difficult to say from remote and without having the full set of modified files (BTW: is is a very bad idea to modify such a class without changing the class name and fix style name as these kinds of changes will break the original command for all normal uses).

On top of this, you have to keep in mind that LAMMPS’ internal interfaces are continually changing due to an ongoing refactoring process (for over two years already, LAMMPS is a large software package) where internal interfaces can change and thus code written for an older version of LAMMPS may not compile anymore. I cannot say whether this applies in this case since there is too little code visible.