Implementation of novel boundary conditions in LAMMPS

Dear Axel and Steve,

It appears that fix property/atom fulfills my purpose. I plan to include the following commands in my input script before defining the simulation box :

fix f1 all property/atom d_rmin d_fsx d_fsy f_fsz

set atom * d_rmin 1000.0 d_fsx 0.0 d_fsy 0.0 d_fsz 0.0

I have modified the void compute of pair_dpd as follows:

//additions

int n = atom->ntypes;

int tmp;

  int index1 = atom->find_custom("rmin", tmp);
  int index2 = atom->find_custom("fsx", tmp);
  int index3 = atom->find_custom("fsy", tmp);
  int index4 = atom->find_custom("fsz", tmp);

  double *rmin = atom->dvector[index1];
  double *fsx = atom->dvector[index2];
  double *fsy = atom->dvector[index3];
  double *fsz = atom->dvector[index4];

  for (i = 0; i < nlocal; i++) {

    rmin[i]=1000.0;
    fsx[i]=0.0;
    fsy[i]=0.0;
    fsz[i]=0.0;
  }

//end of additions

  // 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];
      factor_dpd = special_lj[sbmask(j)];
      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 < 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 = delx*delvx + dely*delvy + 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]*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;

// additions

        if (jtype==n) {

        fsx[i] += delx*fpair;
        fsy[i] += dely*fpair;
        fsz[i] += delz*fpair;
        if (r < rmin[i]) rmin[i] = r ;

        }

// end of additions

        if (newton_pair || j < nlocal) {
          f[j][0] -= delx*fpair;
          f[j][1] -= dely*fpair;
          f[j][2] -= delz*fpair;
//additions
          if (itype==n){
          fsx[j] += delx*fpair;
          fsy[j] += dely*fpair;
          fsz[j] += delz*fpair;
          if (r < rmin[j]) rmin[j] = r ;
// end of additions

I have 3 types of atoms. I expect that the above code will extract the total pairwise DPD force exerted upon each atom by atoms of type 3 and store the results in the vectors fsx, fsy and fsz, which can then be used by my custom fix to modify the velocities at the end of the time-step.

However, some of the previous conversations in this mailing list give me cause for concern. In May 2011, someone asked whether it was possible to output the pairwise force between selected atom types (http://lammps.sandia.gov/threads/msg20801.html) and Axel replied in the negative. As this is precisely what I am trying to do, can I check whether this limitation still exists in the 11Aug17 version of LAMMPS?

A workaround was proposed on this page (https://sjbyrnes.com/OutputElectricForceInLAMMPS.html) but the author mentioned that it was necessary to turn newton off to obtain correct results. Is this necessary in my case? Are the communication routines in fix property/atom sufficient to ensure the accuracy of the force data in fsx, fsy and fsz or do I need to make additions to the pack_reverse and unpack_reverse routines in the atom_vec files?

Best,

Karthik.