git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@15276 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/neigh_derive.cpp

@@ -11,12 +11,23 @@
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
+#include <string.h>
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "atom.h"
+#include "domain.h"
+#include "fix_shear_history.h"
 #include "my_page.h"
 #include "error.h"
 
+
+
+
+// DEBUG
+#include "update.h"
+#include "comm.h"
+
+
 using namespace LAMMPS_NS;
 
 /* ----------------------------------------------------------------------
@@ -222,16 +233,28 @@ void Neighbor::skip_from(NeighList *list)
 /* ----------------------------------------------------------------------
    build skip list for subset of types from parent list
    iskip and ijskip flag which atom types and type pairs to skip
-   this is for granular lists with history, copy the history values from parent
+   if list requests it, preserve shear history via fix shear/history 
 ------------------------------------------------------------------------- */
 
 void Neighbor::skip_from_granular(NeighList *list)
 {
-  int i,j,ii,jj,n,nn,itype,jnum,joriginal;
+  int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes;
+  tagint jtag;
   int *neighptr,*jlist,*touchptr,*touchptr_skip;
   double *shearptr,*shearptr_skip;
 
+  NeighList *listgranhistory;
+  int *npartner;
+  tagint **partner;
+  double (**shearpartner)[3];
+  int **firsttouch;
+  double **firstshear;
+  MyPage<int> *ipage_touch;
+  MyPage<double> *dpage_shear;
+
+  tagint *tag = atom->tag;
   int *type = atom->type;
+  int nlocal = atom->nlocal;
 
   int *ilist = list->ilist;
   int *numneigh = list->numneigh;
@@ -241,23 +264,33 @@ void Neighbor::skip_from_granular(NeighList *list)
   int *ilist_skip = list->listskip->ilist;
   int *numneigh_skip = list->listskip->numneigh;
   int **firstneigh_skip = list->listskip->firstneigh;
-  int **firsttouch_skip = list->listskip->listgranhistory->firstneigh;
-  double **firstshear_skip = list->listskip->listgranhistory->firstdouble;
   int inum_skip = list->listskip->inum;
 
   int *iskip = list->iskip;
   int **ijskip = list->ijskip;
 
-  NeighList *listgranhistory = list->listgranhistory;
-  int **firsttouch = listgranhistory->firstneigh;
-  double **firstshear = listgranhistory->firstdouble;
-  MyPage<int> *ipage_touch = listgranhistory->ipage;
-  MyPage<double> *dpage_shear = listgranhistory->dpage;
+  FixShearHistory *fix_history = list->fix_history;
+  if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nlocal + atom->nghost;
+    npartner = fix_history->npartner;
+    partner = fix_history->partner;
+    shearpartner = fix_history->shearpartner;
+    listgranhistory = list->listgranhistory;
+    firsttouch = listgranhistory->firstneigh;
+    firstshear = listgranhistory->firstdouble;
+    ipage_touch = listgranhistory->ipage;
+    dpage_shear = listgranhistory->dpage;
+    dnum = listgranhistory->dnum;
+    dnumbytes = dnum * sizeof(double);
+  }
 
   int inum = 0;
   ipage->reset();
-  ipage_touch->reset();
-  dpage_shear->reset();
+  if (fix_history) {
+    ipage_touch->reset();
+    dpage_shear->reset();
+  }
 
   // loop over atoms in other list
   // skip I atom entirely if iskip is set for type[I]
@@ -268,15 +301,16 @@ void Neighbor::skip_from_granular(NeighList *list)
     itype = type[i];
     if (iskip[itype]) continue;
 
-    n = nn = 0;
+    n = 0;
     neighptr = ipage->vget();
-    touchptr = ipage_touch->vget();
-    shearptr = dpage_shear->vget();
+    if (fix_history) {
+      nn = 0;
+      touchptr = ipage_touch->vget();
+      shearptr = dpage_shear->vget();
+    }
 
-    // loop over parent non-skip granular list and its history info
+    // loop over parent non-skip granular list and optionally its history info
 
-    touchptr_skip = firsttouch_skip[i];
-    shearptr_skip = firstshear_skip[i];
     jlist = firstneigh_skip[i];
     jnum = numneigh_skip[i];
 
@@ -285,10 +319,166 @@ void Neighbor::skip_from_granular(NeighList *list)
       j = joriginal & NEIGHMASK;
       if (ijskip[itype][type[j]]) continue;
       neighptr[n] = joriginal;
-      touchptr[n++] = touchptr_skip[jj];
-      shearptr[nn++] = shearptr_skip[3*jj];
-      shearptr[nn++] = shearptr_skip[3*jj+1];
-      shearptr[nn++] = shearptr_skip[3*jj+2];
+
+      // no numeric test for current touch
+      // just use FSH partner list to infer it
+      // would require distance calculation for spheres
+      // more complex calculation for surfs
+
+      if (fix_history) {
+        jtag = tag[j];
+        for (m = 0; m < npartner[i]; m++)
+          if (partner[i][m] == jtag) break;
+        if (m < npartner[i]) {
+          touchptr[n] = 1;
+          memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes);
+          nn += dnum;
+        } else {
+          touchptr[n] = 0;
+          memcpy(&shearptr[nn],zeroes,dnumbytes);
+          nn += dnum;
+        }
+      }
+
+      n++;
+    }
+
+    ilist[inum++] = i;
+    firstneigh[i] = neighptr;
+    numneigh[i] = n;
+    ipage->vgot(n);
+    if (ipage->status())
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+
+    if (fix_history) {
+      firsttouch[i] = touchptr;
+      firstshear[i] = shearptr;
+      ipage_touch->vgot(n);
+      dpage_shear->vgot(nn);
+    }
+  }
+
+  list->inum = inum;
+}
+
+/* ----------------------------------------------------------------------
+   build skip list for subset of types from parent list
+   iskip and ijskip flag which atom types and type pairs to skip
+   parent non-skip list used newton off, this skip list is newton on
+   if list requests it, preserve shear history via fix shear/history 
+------------------------------------------------------------------------- */
+
+void Neighbor::skip_from_granular_off2on(NeighList *list)
+{
+  int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,dnum,dnumbytes;
+  tagint itag,jtag;
+  int *neighptr,*jlist,*touchptr,*touchptr_skip;
+  double *shearptr,*shearptr_skip;
+
+  NeighList *listgranhistory;
+  int *npartner;
+  tagint **partner;
+  double (**shearpartner)[3];
+  int **firsttouch;
+  double **firstshear;
+  MyPage<int> *ipage_touch;
+  MyPage<double> *dpage_shear;
+
+  tagint *tag = atom->tag;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  MyPage<int> *ipage = list->ipage;
+
+  int *ilist_skip = list->listskip->ilist;
+  int *numneigh_skip = list->listskip->numneigh;
+  int **firstneigh_skip = list->listskip->firstneigh;
+  int inum_skip = list->listskip->inum;
+
+  int *iskip = list->iskip;
+  int **ijskip = list->ijskip;
+
+  FixShearHistory *fix_history = list->fix_history;
+  if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nlocal + atom->nghost;
+    npartner = fix_history->npartner;
+    partner = fix_history->partner;
+    shearpartner = fix_history->shearpartner;
+    listgranhistory = list->listgranhistory;
+    firsttouch = listgranhistory->firstneigh;
+    firstshear = listgranhistory->firstdouble;
+    ipage_touch = listgranhistory->ipage;
+    dpage_shear = listgranhistory->dpage;
+    dnum = listgranhistory->dnum;
+    dnumbytes = dnum * sizeof(double);
+  }
+
+  int inum = 0;
+  ipage->reset();
+  if (fix_history) {
+    ipage_touch->reset();
+    dpage_shear->reset();
+  }
+
+  // loop over atoms in other list
+  // skip I atom entirely if iskip is set for type[I]
+  // skip I,J pair if ijskip is set for type[I],type[J]
+
+  for (ii = 0; ii < inum_skip; ii++) {
+    i = ilist_skip[ii];
+    itype = type[i];
+    if (iskip[itype]) continue;
+    itag = tag[i];
+
+    n = 0;
+    neighptr = ipage->vget();
+    if (fix_history) {
+      nn = 0;
+      touchptr = ipage_touch->vget();
+      shearptr = dpage_shear->vget();
+    }
+
+    // loop over parent non-skip granular list and optionally its history info
+
+    jlist = firstneigh_skip[i];
+    jnum = numneigh_skip[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      joriginal = jlist[jj];
+      j = joriginal & NEIGHMASK;
+      if (ijskip[itype][type[j]]) continue;
+
+      // only keep I,J when J = ghost if Itag < Jtag
+
+      jtag = tag[j];
+      if (j >= nlocal && jtag < itag) continue;
+
+      neighptr[n] = joriginal;
+
+      // no numeric test for current touch
+      // just use FSH partner list to infer it
+      // would require distance calculation for spheres
+      // more complex calculation for surfs
+
+      if (fix_history) {
+        for (m = 0; m < npartner[i]; m++)
+          if (partner[i][m] == jtag) break;
+        if (m < npartner[i]) {
+          touchptr[n] = 1;
+          memcpy(&shearptr[nn],shearpartner[i][m],dnumbytes);
+          nn += dnum;
+        } else {
+          touchptr[n] = 0;
+          memcpy(&shearptr[nn],zeroes,dnumbytes);
+          nn += dnum;
+        }
+      }
+
+      n++;
     }
 
     ilist[inum++] = i;
@@ -298,10 +488,206 @@ void Neighbor::skip_from_granular(NeighList *list)
     if (ipage->status())
       error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
 
-    firsttouch[i] = touchptr;
-    firstshear[i] = shearptr;
-    ipage_touch->vgot(n);
-    dpage_shear->vgot(nn);
+    if (fix_history) {
+      firsttouch[i] = touchptr;
+      firstshear[i] = shearptr;
+      ipage_touch->vgot(n);
+      dpage_shear->vgot(nn);
+    }
+  }
+
+  list->inum = inum;
+}
+
+/* ----------------------------------------------------------------------
+   build skip list for subset of types from parent list
+   iskip and ijskip flag which atom types and type pairs to skip
+   parent non-skip list used newton off and was not onesided,
+     this skip list is newton on and onesided
+   if list requests it, preserve shear history via fix shear/history 
+------------------------------------------------------------------------- */
+
+void Neighbor::skip_from_granular_off2on_onesided(NeighList *list)
+{
+  int i,j,ii,jj,m,n,nn,itype,jnum,joriginal,flip,dnum,dnumbytes,tmp;
+  tagint itag,jtag;
+  int *surf,*neighptr,*jlist;
+
+  NeighList *listgranhistory;
+  int *npartner;
+  tagint **partner;
+  double (**shearpartner)[3];
+  int **firsttouch;
+  double **firstshear;
+  MyPage<int> *ipage_touch;
+  MyPage<double> *dpage_shear;
+
+  tagint *tag = atom->tag;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+
+  int *ilist = list->ilist;
+  int *numneigh = list->numneigh;
+  int **firstneigh = list->firstneigh;
+  MyPage<int> *ipage = list->ipage;
+
+  int *ilist_skip = list->listskip->ilist;
+  int *numneigh_skip = list->listskip->numneigh;
+  int **firstneigh_skip = list->listskip->firstneigh;
+  int inum_skip = list->listskip->inum;
+
+  int *iskip = list->iskip;
+  int **ijskip = list->ijskip;
+
+  if (domain->dimension == 2) surf = atom->line;
+  else surf = atom->tri;
+
+  FixShearHistory *fix_history = list->fix_history;
+  if (fix_history) {
+    fix_history->nlocal_neigh = nlocal;
+    fix_history->nall_neigh = nlocal + atom->nghost;
+    npartner = fix_history->npartner;
+    partner = fix_history->partner;
+    shearpartner = fix_history->shearpartner;
+    listgranhistory = list->listgranhistory;
+    firsttouch = listgranhistory->firstneigh;
+    firstshear = listgranhistory->firstdouble;
+    ipage_touch = listgranhistory->ipage;
+    dpage_shear = listgranhistory->dpage;
+    dnum = listgranhistory->dnum;
+    dnumbytes = dnum * sizeof(double);
+  }
+
+  int inum = 0;
+  ipage->reset();
+  if (fix_history) {
+    ipage_touch->reset();
+    dpage_shear->reset();
+  }
+
+  // two loops over parent list required, one to count, one to store
+  // because onesided constraint means pair I,J may be stored with I or J
+  // so don't know in advance how much space to alloc for each atom's neighs
+
+  // first loop over atoms in other list to count neighbors
+  // skip I atom entirely if iskip is set for type[I]
+  // skip I,J pair if ijskip is set for type[I],type[J]
+
+  for (i = 0; i < nlocal; i++) numneigh[i] = 0;
+
+  for (ii = 0; ii < inum_skip; ii++) {
+    i = ilist_skip[ii];
+    itype = type[i];
+    if (iskip[itype]) continue;
+    itag = tag[i];
+
+    n = 0;
+
+    // loop over parent non-skip granular list
+
+    jlist = firstneigh_skip[i];
+    jnum = numneigh_skip[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      joriginal = jlist[jj];
+      j = joriginal & NEIGHMASK;
+      if (ijskip[itype][type[j]]) continue;
+
+      // flip I,J if necessary to satisfy onesided constraint
+      // do not keep if I is now ghost
+
+      if (surf[i] >= 0) {
+        if (j >= nlocal) continue;
+        tmp = i;
+        i = j;
+        j = tmp;
+        flip = 1;
+      } else flip = 0;
+
+      numneigh[i]++;
+      if (flip) i = j;
+    }
+  }
+
+  // allocate all per-atom neigh list chunks, including history
+
+  for (i = 0; i < nlocal; i++) {
+    if (numneigh[i] == 0) continue;
+    n = numneigh[i];
+    firstneigh[i] = ipage->get(n);
+    if (ipage->status())
+      error->one(FLERR,"Neighbor list overflow, boost neigh_modify one");
+    if (fix_history) {
+      firsttouch[i] = ipage_touch->get(n);
+      firstshear[i] = dpage_shear->get(dnum*n);
+    }
+  }
+
+  // second loop over atoms in other list to store neighbors
+  // skip I atom entirely if iskip is set for type[I]
+  // skip I,J pair if ijskip is set for type[I],type[J]
+
+  for (i = 0; i < nlocal; i++) numneigh[i] = 0;
+
+  for (ii = 0; ii < inum_skip; ii++) {
+    i = ilist_skip[ii];
+    itype = type[i];
+    if (iskip[itype]) continue;
+    itag = tag[i];
+
+    // loop over parent non-skip granular list and optionally its history info
+
+    jlist = firstneigh_skip[i];
+    jnum = numneigh_skip[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      joriginal = jlist[jj];
+      j = joriginal & NEIGHMASK;
+      if (ijskip[itype][type[j]]) continue;
+
+      // flip I,J if necessary to satisfy onesided constraint
+      // do not keep if I is now ghost
+
+      if (surf[i] >= 0) {
+        if (j >= nlocal) continue;
+        tmp = i;
+        i = j;
+        j = tmp;
+        flip = 1;
+      } else flip = 0;
+
+      // store j in neigh list, not joriginal, like other neigh methods
+      // OK, b/c there is no special list flagging for surfs
+
+      firstneigh[i][numneigh[i]] = j;
+
+      // no numeric test for current touch
+      // just use FSH partner list to infer it
+      // would require complex calculation for surfs
+
+      if (fix_history) {
+        jtag = tag[j];
+        n = numneigh[i];
+        nn = dnum*n;
+        for (m = 0; m < npartner[i]; m++)
+          if (partner[i][m] == jtag) break;
+        if (m < npartner[i]) {
+          firsttouch[i][n] = 1;
+          memcpy(&firstshear[i][nn],shearpartner[i][m],dnumbytes);
+        } else {
+          firsttouch[i][n] = 0;
+          memcpy(&firstshear[i][nn],zeroes,dnumbytes);
+        }
+      }
+
+      numneigh[i]++;
+      if (flip) i = j;
+    }
+
+    // only add atom I to ilist if it has neighbors
+    // fix shear/history allows for this in pre_exchange_onesided()
+
+    if (numneigh[i]) ilist[inum++] = i;
   }
 
   list->inum = inum;

src/neighbor.h

@@ -174,6 +174,8 @@ class Neighbor : protected Pointers {
   int *glist;                  // lists to grow atom arrays every reneigh
   int *slist;                  // lists to grow stencil arrays every reneigh
 
+  double *zeroes;              // vector of zeroes for shear history init
+
   // USER-DPD package
 
   int *bins_ssa;             // ptr to next atom in each bin used by SSA
@@ -236,17 +238,12 @@ class Neighbor : protected Pointers {
   void full_bin_ghost(class NeighList *);
   void full_multi(class NeighList *);
 
-  void half_from_full_no_newton(class NeighList *);
-  void half_from_full_newton(class NeighList *);
-  void skip_from(class NeighList *);
-  void skip_from_granular(class NeighList *);
-  void skip_from_respa(class NeighList *);
-  void copy_from(class NeighList *);
-
   void granular_nsq_no_newton(class NeighList *);
   void granular_nsq_newton(class NeighList *);
+  void granular_nsq_newton_onesided(class NeighList *);
   void granular_bin_no_newton(class NeighList *);
   void granular_bin_newton(class NeighList *);
+  void granular_bin_newton_onesided(class NeighList *);
   void granular_bin_newton_tri(class NeighList *);
 
   void respa_nsq_no_newton(class NeighList *);
@@ -255,6 +252,15 @@ class Neighbor : protected Pointers {
   void respa_bin_newton(class NeighList *);
   void respa_bin_newton_tri(class NeighList *);
 
+  void half_from_full_no_newton(class NeighList *);
+  void half_from_full_newton(class NeighList *);
+  void skip_from(class NeighList *);
+  void skip_from_granular(class NeighList *);
+  void skip_from_granular_off2on(class NeighList *);
+  void skip_from_granular_off2on_onesided(class NeighList *);
+  void skip_from_respa(class NeighList *);
+  void copy_from(class NeighList *);
+
   // include prototypes for multi-threaded neighbor lists
   // builds or their corresponding dummy versions
 

src/pair.h

@@ -170,7 +170,7 @@ class Pair : protected Pointers {
   virtual void write_data(FILE *) {}
   virtual void write_data_all(FILE *) {}
 
-  virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;}
+  virtual int pack_forward_comm(int, int *, double *, int, int *) {printf("WHAT\n"); return 0;}
   virtual void unpack_forward_comm(int, int, double *) {}
   virtual int pack_forward_comm_kokkos(int, DAT::tdual_int_2d, 
                                        int, DAT::tdual_xfloat_1d &, 

src/pair_hybrid.cpp

@@ -491,7 +491,7 @@ void PairHybrid::init_style()
 
   for (istyle = 0; istyle < nstyles; istyle++) styles[istyle]->init_style();
 
-  // create skip lists for each pair neigh request
+  // create skip lists inside each pair neigh request
   // any kind of list can have its skip flag set at this stage
 
   for (i = 0; i < neighbor->nrequest; i++) {
@@ -539,6 +539,7 @@ void PairHybrid::init_style()
     // if any skipping occurs
     // set request->skip and copy iskip and ijskip into request
     // else delete iskip and ijskip
+    // no skipping if pair style assigned to all type pairs
 
     skip = 0;
     for (itype = 1; itype <= ntypes; itype++)
@@ -618,7 +619,8 @@ void PairHybrid::setup()
 }
 
 /* ----------------------------------------------------------------------
-   combine sub-style neigh list requests and create new ones if needed
+   examine sub-style neigh list requests
+   create new parent requests if needed, to derive sub-style requests from
 ------------------------------------------------------------------------- */
 
 void PairHybrid::modify_requests()
@@ -626,29 +628,37 @@ void PairHybrid::modify_requests()
   int i,j;
   NeighRequest *irq,*jrq;
 
-  // loop over pair requests only
-  // if list is skip list and not copy, look for non-skip list of same kind
-  // if one exists, point at that one via otherlist
-  // else make new non-skip request of same kind and point at that one
-  //   don't bother to set ID for new request, since pair hybrid ignores list
-  // only exception is half_from_full:
-  //   ignore it, turn off skip, since it will derive from its skip parent
-  // after possible new request creation, unset skip flag and otherlist
-  //   for these derived lists: granhistory, rRESPA inner/middle
-  //   this prevents neighbor from treating them as skip lists
-  // copy list check is for pair style = hybrid/overlay
-  //   which invokes this routine
+  // loop over pair requests only, including those added during looping
+
+  int nrequest_original = neighbor->nrequest;
 
   for (i = 0; i < neighbor->nrequest; i++) {
     if (!neighbor->requests[i]->pair) continue;
 
+    // nothing more to do if this request:
+    //   is not a skip list
+    //   is a copy or half_from_full or granhistory list
+    // copy list setup is from pair style = hybrid/overlay
+    //   which invokes this method at end of its modify_requests()
+    // if granhistory, turn off skip, since each gran sub-style
+    //   its own history list, parent gran list does not have history
+    // if half_from_full, turn off skip, since it will derive 
+    //   from its full parent and its skip status
+
     irq = neighbor->requests[i];
-    if (irq->skip == 0 || irq->copy) continue;
-    if (irq->half_from_full) {
+    if (irq->skip == 0) continue;
+    if (irq->copy) continue;
+    if (irq->granhistory || irq->half_from_full) {
       irq->skip = 0;
       continue;
     }
 
+    // look for another list that matches via same_kind() and is not a skip list
+    // if one exists, point at that one via otherlist
+    // else make new parent request via copy_request() and point at that one
+    //   new parent list is not a skip list
+    //   parent does not need its ID set, since pair hybrid does not use it
+
     for (j = 0; j < neighbor->nrequest; j++) {
       if (!neighbor->requests[j]->pair) continue;
       jrq = neighbor->requests[j];
@@ -662,11 +672,44 @@ void PairHybrid::modify_requests()
       irq->otherlist = newrequest;
     }
 
-    if (irq->granhistory || irq->respainner || irq->respamiddle) {
+    // for rRESPA inner/middle lists,
+    //   which just created or set otherlist to parent:
+    // unset skip flag and otherlist
+    //   this prevents neighbor from treating them as skip lists
+
+    if (irq->respainner || irq->respamiddle) {
       irq->skip = 0;
       irq->otherlist = -1;
     }
   }
+
+  // adjustments to newly added granular parent requests (gran = 1)
+  // parent newton = 2 if has children with granonesided = 0 and 1
+  //   else newton = 0 = setting of children
+  // parent gran onesided = 0 if has children with granonesided = 0 and 1
+  //   else onesided = setting of children
+
+  for (i = nrequest_original; i < neighbor->nrequest; i++) {
+    if (!neighbor->requests[i]->pair) continue;
+    if (!neighbor->requests[i]->gran) continue;
+    irq = neighbor->requests[i];
+
+    int onesided = -1;
+    for (j = 0; j < nrequest_original; j++) {
+      if (!neighbor->requests[j]->pair) continue;
+      if (!neighbor->requests[j]->gran) continue;
+      jrq = neighbor->requests[j];
+      
+      if (onesided < 0) onesided = jrq->granonesided;
+      else if (onesided != jrq->granonesided) onesided = 2;
+      if (onesided == 2) break;
+    }
+
+    if (onesided == 2) {
+      irq->newton = 2;
+      irq->granonesided = 0;
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------